THE  USE  OF  UNDETERMINED  MULTIPLIERS  IN  THE 
ANALYSIS  AND  OPTIMUM  DESIGN  OF  ELASTIC  STRUCTURES 


BY 

TZU-YANG  CHEN 


A DISSERTATION  PRESENTED  TO  THE  GRADUATE  SCHOOL 
OF  THE  UNIVERSITY  OF  FLORIDA  IN  PARTIAL  FULFILLMENT 
OF  THE  REQUIREMENTS  FOR  THE  DEGREE  OF 
DOCTOR  OF  PHILOSOPHY 

UNIVERSITY  OF  FLORIDA 


1989 


ACKNOWLEDGEMENTS 


The  author  sincerely  wishes  to  thank  the  chairman  of  his 
supervisory  committee,  Dr.  Ali  A.  Seireg,  for  his  guidance  and 
advice  throughout  the  research. 

The  author  wishes  to  express  his  appreciation  to  the 
other  members  of  his  supervisory  committee.  Dr.  Joseph  Duffy, 
Dr.  George  N.  Sandor,  Dr.  Gary  K.  Matthew  and  Dr.  Ulrich  H. 
Kurzweg  for  their  assistance  in  this  work. 

The  author  would  like  to  thank  Wei  Lin  and  Chun-Yin  Wu 
for  their  help  in  preparing  the  drawings  and  assistance  in 
using  computers,  and  Lindsay  Wells  for  his  help  in  getting 
this  dissertation  in  proper  form.  Without  their  help,  this 
work  could  not  be  done. 

The  author  also  wants  to  thank  his  friends  for  helping 
him  to  get  accustomed  to  life  in  the  United  States  of  America. 

Finally,  the  author  would  like  to  express  his  deepest 
appreciation  to  his  parents  and  his  brother  and  sister  for 
their  financial  support,  encouragement  and  enduring  patience 
throughout  this  study. 


ii 


TABLE  OF  CONTENTS 


ACKNOWLEDGEMENTS  ii 

ABSTRACT  V 

CHAPTER 

1 INTRODUCTION  1 

The  Energy  Principle  1 

Rayleigh-Ritz  Method  Using  Undetermined 

Multipliers  in  Structural  Mechanics  4 

2 THEORETICAL  BACKGROUND 8 

Principle  of  Virtual  Work  8 

Principle  of  Minimum  Potential  Energy  10 

Application  of  Variational  Calculus 

in  Structural  Mechanics  12 

Rayleigh-Ritz  Method  13 

Lagrangian  Multiplier  Method  18 

Rayleigh-Ritz  Method  Using 

Lagrangian  Multipliers  21 

3 APPLICATIONS  TO  BEAM  PROBLEMS 23 

Introduction  23 

Simply-Supported  Beams  25 

Clamped-Supported  Beams  27 

Clamped-Clamped  Beams  29 

Cantilever  Beams  32 

Cantilever  Beams  with  the 

Free  End  Flexibly  Supported 3 6 

Beams  Supported  on  Elastic  Foundations 41 

Multi-Span  Beams  42 

4 APPLICATIONS  TO  PLATE  PROBLEMS  47 

Introduction  47 

Plates  with  All  Edges  Simply  Supported  50 

Plates  with  One  Edge  Clamped 

and  Three  Edges  Simply  Supported  55 

Plates  with  All  Edges  Clamped 62 

Nonrectangular  Plates  66 

iii 


Plates  with  Mixed  Boundary  Conditions  70 

Stiffened  Plates  74 

Honeycomb  Sandwich  Plates  80 

5 APPLICATIONS  TO  PLATES  WITH  HOLES  90 

Introduction  90 

Square  Plates  with  Clamped  Square  Holes  92 

Simulation  of  a Cantilever  Beam 

From  a Clamped  Beam 93 

Square  Plates  with  Free  Square  Holes  101 

6 APPLICATIONS  TO  SHELL  PROBLEMS  107 

Introduction  107 

Cylindrical  Shells  of  Infinite  Length 

Supported  by  a Rigid  Ring 109 

Cylindrical  Shells  of  Infinite  Length 

Reinforced  by  a Narrow  Ring  115 

Cylindrical  Shells  of  Finite  Length 

with  Both  Ends  Clamped 121 

Clamped  Cylindrical  Shells  with 

Circumferential  Ring  Reinforcements  125 

Shells  Reinforced  by  Rings 

and  Longitudinal  Stiffeners  127 

7 APPLICATIONS  TO  STRUCTURAL  DESIGN  136 

Beam  Design 137 

Plate  Design 148 

8 CONCLUSIONS  AND  RECOMMENDATIONS  165 

Conclusions 165 

Selection  of  the  Coordinate  Functions  167 

Limitations  of  the  Method 168 

Recommendations  169 

APPENDIX  170 

REFERENCES 192 

BIOGRAPHICAL  SKETCH  196 


iv 


Abstract  of  Dissertation  Presented  to  the  Graduate  School 
of  the  University  of  Florida  in  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of  Doctor  of  Philosophy 

THE  USE  OF  UNDETERMINED  MULTIPLIERS  IN  THE 
ANALYSIS  AND  OPTIMUM  DESIGN  OF  ELASTIC  STRUCTURES 

By 

Tzu-Yang  Chen 
August  1989 

Chairman:  Dr.  Ali  A.  Seireg 

Major  Department:  Mechanical  Engineering 

The  study  reported  in  this  dissertation  investigates  the 

use  of  undetermined  multipliers  with  the  Rayleigh-Ritz  method 

for  the  solution  of  general  elasticity  problems  with  complex 

geometries  and  boundary  conditions.  Beam,  plate  and  shell 

problems  are  considered  for  illustrations.  A procedure 

incorporating  the  undetermined  multiplier  method  is  given  for 

the  optimum  design  of  beams  and  stiffened  plates. 

Several  cases  of  beams  on  multiple  supports,  plates  with 

holes,  plates  with  reinforcements  and  reinforced  shells  are 

investigated  through  the  use  of  the  Raleigh-Ritz  method 

incorporating  the  undetermined  multipliers.  The  result  shows 

that  this  approach  provides  an  excellent  tool  for  the  analysis 

and  optimum  design  of  complex  structures. 

Although  the  method  requires  some  manual  integration,  the 

result  of  the  integration  can  be  applied  to  similar  problems 


V 


with  different  boundary  conditions  without  any  major  change 
in  formulation.  The  formulation  of  this  method  has  a unified 
character  for  similar  problems.  The  simplicity  and  elegance 
of  this  method  is  particularly  evident  when  dealing  with  a 
class  of  problems  with  different  boundary  conditions. 

This  method  gives  a positive  definite  matrix  for  solving 
the  generated  system  of  linear  equations.  Choleski's 
decomposition  scheme  is  incorporated  to  further  increase  the 
overall  efficiency  of  the  solution.  The  Newton-Raphson  method 
is  also  incorporated  with  the  undetermined  multiplier  approach 
to  give  an  excellent  and  efficient  design  procedure.  This 
method  shows  considerable  advantages  over  the  finite  element 
method  for  the  class  of  problems  considered  in  this  study. 
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CHAPTER  1 
INTRODUCTION 

The  Energy  Principle 

Structural  mechanics  has  long  been  a well-developed 
branch  of  science.  Many  principles  have  been  formulated  to 
describe  the  stresses  and  deformations  in  elastic  bodies.  In 
applying  these  principles  to  the  analysis  and  design  of 
mechanical  system  components,  an  essential  step  is  to 
mathematically  model  the  real  situation,  and  obtain  either  an 
analytical  or  numerical  solution  from  the  model. 

From  the  point  of  view  of  formulating  theoretical  models 
for  the  analysis  of  elastic  structures,  and  the  determination 
of  numerical  solutions,  the  energy  principle  holds  a position 
of  vital  importance[l] [2] . The  essence  of  the  energy 
principle  for  mechanical  structures  is  that  when  the  behavior 
of  an  elastic  element  follows  the  rule  of  "minimum  energy," 
the  application  of  variational  calculus  will  result  in  a 
differential  eguation  governing  the  behavior  of  the  element. 

There  are  many  methods  which  have  been  developed  on  the 
basis  of  the  energy  principle.  Most  of  these  methods  are 
based  on  the  "principle  of  minimum  potential  energy."  The 
energy  principle  can  be  applied  by  using  the  formal  procedure 
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of  the  calculus  of  variations  to  determine  the  stationary 
condition  for  a functional  (the  potential  energy  in  the  case 
of  elastic  systems)  so  as  to  yield  the  differential  equations 
and  appropriate  boundary  conditions [3] [4] [5] [6] . Instead  of 
using  the  formal  procedure,  however,  there  are  alternatives. 
These  are  the  direct  methods [1],  which  substitute  an  assumed 
solution  involving  adjustable  parameters  for  the  functional. 
The  stationary  condition  is  determined  with  respect  to  these 
parameters.  Among  these  methods,  the  most  widely  used  are  the 
Rayleigh-Ritz  method[l] [7] , the  Galerkin  method [1],  and  the 
finite  element  method[8] [9] , which  is  the  most  popular 
approach.  All  of  these  methods  have  inherent  advantages  and 
disadvantages,  depending  upon  the  particular  application. 

The  finite  element  method  is  very  flexible [ 10 ] . It  can 
be  applied  to  general  elasticity  problems,  even  in  nonlinear 
domains [ 11] [12 ] . It  can  even  be  applied  to  plasticity 
problems [ 13 ] . This  method,  however,  has  disadvantages,  the 
most  severe  being  the  extensive  preprocessing  and 
postprocessing  effort  that  is  required.  Another  disadvantage 
is  that  it  is  necessary  to  choose  a very  fine  mesh  in  regions 
with  steep  stress  gradients [7 ] . This  leads  to  the 
construction  and  solution  of  a large  number  of  equations  with 
a large  number  of  unknowns.  In  other  words,  if  there  is  no 
efficient  computing  machine,  the  finite  element  method  cannot 
be  practically  used.  There  are  also  certain  difficulties  when 
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the  finite  element  method  is  applied  to  problems  including 
curved  boundaries [7 ] . 

Basically,  with  all  boundary  conditions  satisfied  in  an 
elasticity  problem,  the  Galerkin  method  directly  provides  a 
solution  to  the  differential  equations,  and  is  applicable 
whether  the  transformation  into  a variational  problem  is 
possible  or  not.  In  other  words,  the  Galerkin  method  can 
supply  an  approximate  solution  without  the  construction  of  the 
functional  (the  potential  energy) . This  advantage,  however, 
is  not  important  in  view  of  the  fact  that  in  elasticity,  an 
expression  for  the  total  potential  energy  can  always  be 
formulated.  Also,  its  limitations,  which  have  been  widely 
discussed  by  Kanagasundaram  and  Harimhallo [ 14 ] , make  it  more 
difficult  to  apply  than  the  Rayleigh-Ritz  method. 

In  contrast  to  the  Galerkin  method,  the  construction  of 
the  functional  is  essential  to  the  Rayleigh-Ritz  method. 
This,  however,  presents  no  difficulty  in  applying  the  method, 
because  the  potential  energy  can  be  formed  directly  by 
differentiation  under  the  integral  sign.  Furthermore,  while 
the  Galerkin  method  requires  the  satisfaction  of  any  traction 
boundary  condition  together  with  the  geometric  boundary 
conditions,  the  Rayleigh-Ritz  method  requires  only  the 
satisfaction  of  the  geometrical  boundary  conditions [3 ] . It 
has  been  pointed  out[l]  that  the  Rayleigh-Ritz  method  is  very 
effective  in  obtaining  approximate  solutions  to  boundary  value 
problems.  Further,  this  method  is  more  efficient  than  the 
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finite  element  method,  since  it  does  not  require  the  solution 
of  as  many  equations.  Actually,  the  Rayleigh-Ritz  method  is 
regarded  as  the  basis  of  the  direct  methods [7],  The  finite 
element  method  can  be  considered  a special  case  of  the 
Rayleigh-Ritz  method. 

Ravleiqh-Ritz  Method  Using  Undetermined 
Multipliers  in  Structural  Mechanics 

The  application  of  the  Rayleigh-Ritz  method  in  simple 
structural  analysis  can  be  found  in  many  books,  for  example 
references [ 3 ] [4]  [5]  [15]  [16]  . These  applications  are  generally 
limited  to  the  simplest  problems,  such  as  a simply-supported 
beam,  or  a simply-supported  plate,  and  only  occasionally  to 
more  complex  problems [17] . This  is  because  the  Rayleigh-Ritz 
method  requires  that  the  assumed  solution  satisfies  the 
boundary  conditions,  and  it  may  not  be  always  possible  to  find 
a solution  which  meets  the  boundary  condition  requirement  and 
minimizes  the  functional  at  the  same  time. 

In  this  study,  an  alternative  way  to  apply  the  Rayleigh- 
Ritz  method  will  be  developed.  The  alternative,  which  will 
be  called  the  modified  Rayleigh-Ritz  method,  incorporates  the 
Lagrangian  rule [18]  [19]  [20]  in  a way  such  that  the 
satisfaction  of  the  boundary  conditions  by  the  assumed 
solution  is  no  longer  necessary.  Through  the  use  of 
Lagrangian  multipliers  (also  known  as  undetermined 
multipliers) , the  boundary  conditions  can  be  forced  to  be 
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satisfied.  The  result  will  still  minimize  the  potential 
energy,  only  this  time  the  minimization  is  done  under  the 
restriction  that  the  boundary  conditions  be  satisfied.  In 
this  way,  it  is  not  necessary  to  find  an  assumed  solution 
which  satisfies  the  boundary  conditions  and  the  minimization 
requirement  of  the  potential  energy  at  the  same  time. 

By  applying  the  modified  Rayleigh-Ritz  method,  many 
difficult  problems  in  classic  elasticity  theory  can  be  solved. 
Several  authors [ 3 ] [21 ] [22 ] have  applied  the  concept  of  the 
Lagrangian  rule  to  handle  constraints  inside  the  domains  of 
plates.  These  applications,  however,  still  require  that  the 
assumed  solution  satisfies  both  the  boundary  conditions  and 
the  minimization  of  the  potential  energy  simultaneously.  It 
is  proposed  that  boundary  conditions,  as  well  as  other 
constraints,  can  be  satisfied  by  modifying  the  Rayleigh-Ritz 
method  with  Lagrangian  multipliers,  and  it  would  not  be 
necessary  to  explicitly  find  an  assumed  solution  which 
satisfies  both  the  boundary  conditions  and  the  minimization 
of  the  potential  energy. 

Although  classical  theories  have  been  very  successful  in 
the  analysis  of  many  structural  elements,  including  beams, 
plates  and  shells,  there  are  still  many  problems  in  these 
categories  that  are  difficult  to  solve  by  classical  methods. 
An  example  is  the  problem  of  a plate  with  mixed  boundaries. 
There  are  generally  no  closed-form  solutions  for  these 
problems,  and  their  analysis  is  left  to  sophisticated  and 
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expensive  numerical  schemes  such  as  the  finite  element  method. 
However,  the  Rayleigh-Ritz  method  incorporating  the 
undetermined  multipliers  can  be  applied  to  solve  this  type  of 
problem.  It  has  been  found  that  this  modified  Rayleigh-Ritz 
method  has  good  potential  for  dealing  with  such  problems  more 
efficiently  than  the  finite  element  method. 

Another  problem  that  the  classical  method  cannot  handle 
easily  is  the  problem  of  a stiffened  plate.  There  are  two 
common  ways  to  deal  with  this  type  of  problem.  One  is  to 
divide  the  plate  into  smaller  elements  bounded  by  the 
stiffeners,  which  generally  reguires  the  solution  of  a large 
system  of  simultaneous  eguations[3] . The  other  approach  is 
to  transfer  the  stiffened  plate  into  an  equivalent  orthotropic 
plate.  This  approach  has  been  favored  by  many 
researchers[23] [24] [25] [26]  for  its  simplicity.  However, 
there  is  no  guarantee  that  the  stresses  obtained  will  match 
at  each  point  of  the  stiffened  plate.  Although  this 
difficulty  can  be  overcome  by  choosing  the  spacing  of  the 
stiffeners  to  be  much  smaller  than  the  dimensions  of  the 
plate [27],  it  does  not  provide  a practical  way  to  handle  more 
general  problems  of  this  type.  On  the  other  hand,  by 
establishing  suitable  constraints  between  the  plate  and  the 
stiffeners,  this  type  of  problem  becomes  well  suited  for 
efficient  and  effective  solution  by  the  modified  Rayleigh-Ritz 


method . 
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The  modified  Rayleigh-Ritz  method  can  also  provide  an 
excellent  tool  for  optimizing  elastic  structures.  In 
structural  optimization,  numerous  technigues  are  used[28], 
such  as  design  sensitivity,  integer  programming,  and 
seguential  linear  programming.  Usually,  the  basic  method  used 
for  the  optimization  of  structures  of  complex  shapes  is  the 
finite  element  method[29].  As  mentioned  before,  the  finite 
element  method  requires  considerable  preprocessing  and 
postprocessing  effort  in  analyzing  structures  and,  thus,  is 
an  expensive  method.  For  example,  for  the  aircraft  wing 
structures  discussed  by  Seong[29],  it  may  take  the  finite 
element  method  many  hours  to  complete  a design  on  a super 
computer.  Furtheirmore,  it  may  be  difficult  to  ensure  the 
accuracy  of  the  finite  element  method  in  some  cases[30]. 

Since  the  optimal  design  process  is  iterative  in  nature 
and  requires  many  evaluations  of  the  objective  function,  the 
overall  efficiency  of  the  proposed  method,  when  used  in 
conjunction  with  efficient  schemes  for  solving  nonlinear 
equations,  makes  it  a valuable  tool  to  use  in  a wide  variety 
of  structural  design  tasks. 


CHAPTER  2 

THEORETICAL  BACKGROUND 
Principle  of  Virtual  Work 

As  mentioned  in  the  previous  chapter,  the  energy 
principle  is  very  important  in  solving  problems  of  elastic 
structures.  It  originates  from  the  principle  of  virtual  work, 
which  states  that  for  a body  in  static  equilibrium,  the 
virtual  work  done  by  the  loads  acting  upon  it  is 
zero[l] [15] [31] . 

Virtual  work  is  defined  as  the  work  done  by  the  loads 
when  the  body  undergoes  a fictitious  displacement  called 
virtual  displacement.  The  virtual  displacement  may  consist 
of  a translation  in  any  direction,  a rotation  about  any  axis, 
or  a combination  of  both,  provided  that  the  translation  and 
the  rotation  are  admissible.  In  mathematic  form,  the 
principle  of  virtual  work  is 

<5W  = 0 (2.1) 

The  symbol  W denotes  work,  and  5W  is  the  virtual  work. 

For  a deformable  body  as  shown  in  figure  2.1,  where  the 
bold  line  represents  the  undeformed  configuration  and  the  thin 
line  represents  the  deformed  configuration,  equation  (2.1)  can 
be  written  as [15] 
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5W  = 6W3,,  + = 0 


(2.2) 


Figure  2.1  A deformable  body 


In  equation  (2.2),  represents  the  work  done  by  the 
external  loads  during  the  process  of  an  admissible  virtual 
deformation  of  the  body.  An  admissible  deformation  is  one 
that  satisfies  the  boundary  conditions.  If  the  body  shown  in 
figure  2.1  is  regarded  as  a combination  of  an  infinite  number 
of  small  elements,  is  the  work  done  by  the  internal 
forces  between  small  elements  of  the  body.  If  Hooke's  law  or 
another  load-deformation  law  is  applicable  to  the  deformable 
body,  the  internal  forces  will  have  a potential  energy  U which 
is  equal  to  the  internal  work  with  opposite  sign[l]: 

U = (2.3) 
Substituting  equation  (2.3)  in  equation  (2.2)  yields 

«5Wext  = (2.4) 
This  means  that  the  work  done  by  external  loads  is  transferred 
to  the  strain  energy  stored  in  the  deformed  body. 
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Principle  of  Minimum  Potential  Energy 

If  the  external  loading  is  conservative,  it  has  a 
potential  ^<[1]: 

^ = -Wext  (2.5) 
Substituting  equations  (2.3)  and  (2.5)  in  equation  (2.2) 
yields 

6W  = ( -5U  ) + ( -<S$  ) = 0 (2.6) 
or  equivalently, 

(SU  + 5$  = 5(  U + ^>  ) =0  (2.7) 
The  sum  of  the  potential  of  both  external  and  internal  forces 
is  defined  as  the  total  potential  of  the  system  and  denoted 
by  n , i . e . , 

n = U + (2.8) 
Then,  from  equation  (2.7), 

(Sn  = 6(  U + $ ) = 0 (2.9) 
This  equation  states  that  for  any  admissible  virtual 
deformation  the  deformable  body  can  take,  the  change  in  total 
potential  energy  is  zero [16]. 

If  the  virtual  deformation  is  regarded  as  an 
infinitesimal  variation  of  the  body's  configuration  from  its 
current  equilibrium  state,  then  611  represents  the  variation 
of  total  potential  energy  due  to  the  change  in  configuration. 
Because  of  equation  (2.9),  it  is  clear  that  the  variation  in 
potential  energy  should  be  zero.  In  other  words,  n has  a 
stationary  value  when  the  body  is  in  equilibrium[ 15] . 
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Equation  (2.9)  can  then  be  interpreted  to  mean  that  of  all 
admissible  configurations  a deformable  body  can  take,  those 
satisfying  the  equilibrium  condition  result  in  a stationary 
value  of  the  total  potential  energy [16].  This  constitutes  the 
principle  of  stationary  potential  energy,  and  an  admissible 
configuration  means  a configuration  satisfying  the  boundary 
conditions . 

By  combining  equation  (2.5)  with  equations  (2.6)  and 
(2.9),  the  principle  of  stationary  potential  energy  can  be 
written  as 

6n  = 5U  - = 5(  U - ) = 0 (2.10) 
If  the  unloaded  configuration  of  the  deformable  body  (the  bold 
line  in  figure  2.1)  is  taken  as  the  reference  datum,  the 
potential  energy  contributed  by  external  loads  is  the  work 
done  by  these  loads  when  the  body  undergoes  deformation  from 
its  unloaded  state  to  the  final  equilibrium  configuration  (the 
thin  line  in  figure  2.1) [15]. 

A stationary  value,  by  definition,  may  be  a minimum,  a 
maximum  or  a neutral  value.  However,  since  the  potential 
energy  function  can  be  shown  to  be  positive  definite[16] , it 
can  be  concluded  that  II  is  a minimum  when  the  body  is  in 
static  equilibrium.  It  can  then  be  stated  that  of  all 
admissible  configurations,  those  satisfying  the  condition  of 
equilibrium  result  in  a minimum  value  of  the  total  potential 
energy  of  the  system [16]. 


Application  of  Variational  Calculus 
in  Structural  Mechanics 
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Variational  calculus  is  mainly  concerned  with  the 
stationary  value  of  a definite  integral  called  a functional. 
The  integral  involves  some  functions  which  are  not  known,  and 
the  purpose  is  to  determine  the  unknown  functions  so  as  to 
make  the  functional  stationary.  In  view  of  the  fact  that  the 
configuration  of  an  elastic  body  in  static  equilibrium  results 
in  a minimum  total  potential  energy  of  the  body,  it  can  be 
understood  that  this  type  of  problem  is  well  suited  to 
variational  calculus.  As  a matter  of  fact,  this  is  a major 
application  of  variational  calculus. 

Consider  the  following  functional: 


n = 


^xl 


F(x,v,v',v")  dx 


(2.11) 


xO 


where  F is  an  explicit  function  of  x,  v,  v',  and  v",  and  Xg 
and  Xj  are  fixed  numbers.  A function  v=v(x)  is  sought  to  make 
n stationary.  If  II  is  an  extremum,  then  511=0.  The  optimum 
function  v which  makes  the  functional  an  extremum  satisfies 
Euler's  equation,  which,  in  this  case,  is [4] [5] 


3F  d dF  d' 

- — (^)  + 


dF 

— T ( 

dx^  3v" 


) = 0 


(2.12) 


3v  dx  ' 3v' 
provided  that  v and  v'  are  prescribed  at  x=Xg  and  x=Xj. 

If  V and  v'  are  not  specified,  the  following  conditions. 


in  addition  to  equation  (2.12),  should  also  be  satisfied  to 
make  511  vanish: 
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and 


5v'  dx  ^ 3v"  ^ 


(2.13) 


3F 

= 0 

3v" 


(2.14) 


Equations  (2.13)  and  (2.14)  are  customarily  referred  to  as 
natural  boundary  conditions [4] [5] . 

Consider  the  problem  of  a beam  subjected  to  a lateral 
load  p(x).  If  the  origin  of  a coordinate  system  is  located 
at  one  end  of  the  beam,  and  the  positive  x axis  of  the 
coordinate  system  is  running  toward  the  other  end,  then  the 
potential  energy  is [15] 


n 


p(x)v]  dx 


(2.15) 


where  f is  the  length  of  the  beam,  El  is  the  bending  rigidity. 
In  equation  (2.15),  v denotes  the  deflection  function  (the 
equilibrium  configuration)  of  the  beam.  Substituting 
(l/2)EI(v")^-pv  as  F in  equation  (2.12)  yields 
d^  dV 

■^(EI-^)  - p(x)  = 0 (2.16) 

which  is  the  well-known  governing  equation  for  a beam 
subjected  to  a lateral  load  function  p(x) . 


Ravleiqh-Ritz  Method 

In  the  previous  sections,  as  an  application  of 
variational  calculus  to  beam  problems,  the  equivalence  between 
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the  differential  ecpiation  which  governs  the  deflection  of  a 
beam  and  the  energy  formulation  was  demonstrated.  This 
eguivalence  can  be  utilized  in  many  ways  to  numerically  handle 
problems  of  elastic  systems.  One  of  these  methods  is  the 
Rayleigh-Ritz  method[l] [5] . 

A deformable  body  with  its  properties  (mass,  stiffness, 
etc.)  distributed  in  a continuous  manner  is  referred  to  as  a 
continuum.  A continuum  usually  has  an  infinite  number  of 
degrees  of  freedom  and,  thus,  can  assume  an  infinite  number 
of  deformed  configurations.  It  is  customary  to  describe  the 
response  of  such  a system  with  a differential  eguation,  or 
several  equations,  together  with  boundary  conditions. 
Sometimes,  it  is  hard  to  solve  the  differential  equation (s) 
with  the  prescribed  boundary  conditions.  It  is,  however, 
possible  to  approximate  the  real  system  with  a system 
possessing  only  a finite  number  of  degrees  of  freedom  by 
assuming  certain  characteristics  of  the  deformation. 

To  do  so,  a reasonable  deformed  configuration  involving 
some  undetermined  parameters  is  assumed.  Taking  a one- 
dimensional problem  as  an  example,  if  v(x)  is  the  function 
which  makes  the  potential  energy  II  minimum,  then  it  is 
approximated  in  the  Rayleigh-Ritz  method  by  another  function 
Vn(x).  This  function  is  a linear  combination  of  a set  of 
prescribed  coordinate  functions  §p(x)  in  the  following  form: 

v(x)  « v,(x)  =1^  a,#Jx) 


(2.17) 
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where  are  unknown  parameters,  called  displacement 

parameters,  and  are  to  be  determined  by  minimizing  the 
potential  energy  II.  Usually,  the  coordinate  functions  are 

chosen  in  such  a way  that  no  matter  how  the  parameters  a^  are 
determined,  expression  (2.17)  always  satisfies  the  boundary 
conditions. 

With  equation  (2.17),  the  potential  energy  II  can  be 
calculated  on  the  basis  of  the  knowledge  of  the  coordinate 
functions  The  potential  energy  becomes  a function  of  the 

displacement  parameters  a^^  only.  Instead  of  dealing  with  a 
functional  as  described  in  the  previous  sections,  we  are  now 
handling  a function  and,  thus,  we  have  an  ordinary  extremum 
problem.  The  potential  energy  is  put  in  the  form  of  a 
function  of  the  parameters  a^,  i.e., 

n n ( a^^ , 3.2  / • • • / (2.18) 

where  the  domain  of  the  function  II  may  be  regarded  as  the 
whole  N-dimensional  Euclidean  space  with  the  parameters  a^  as 
the  set  of  coordinates.  The  potential  energy  function  II  is, 
at  least,  of  the  class  in  the  N-dimensional  Euclidean 
space[2].  This  means  that  the  function  itself,  as  well  as 
its  first  derivative,  is  continuous. 

In  accordance  with  differential  calculus,  the  necessary 
and  sufficient  condition  for  a stationary  point  is 

an 

T = 0 n=l,2,...,N  (2.19) 

da 

The  solution  of  the  simultaneous  algebraic  equations  (2.19) 
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corresponds  to  the  stationary  point  of  the  function  II.  It 
has  already  been  shown  that  this  stationary  point  is  a minimum 
point. 

In  linear  elasticity  problems,  II  is  quadratic  in  a^  so 
the  differentiation  indicated  (2.19)  produces  a set  of  linear 
simultaneous  equations [1] . This  considerably  reduces  the  work 
to  solve  the  problem.  However,  it  must  be  understood  that  the 
variational  problem  has  been  replaced  by  an  ordinary  extremum 
problem  through  the  use  of  the  set  of  prescribed  coordinate 
functions  and,  in  general,  these  two  are  not  equivalent. 
Therefore,  the  solution  obtained  from  the  ordinary  extremum 
problem  will  only  be  an  approximation  to  that  of  the 
variational  problem. 

The  efficiency  of  the  Rayleigh-Ritz  method  is  clearly 
dependent  upon  the  choice  of  the  coordinate  functions  §^(x). 
In  the  sense  of  efficiency,  it  is  advantageous  to  satisfy  all 
the  boundary  conditions  in  applying  this  procedure.  This  is, 
however,  not  necessary  in  the  sense  of  the  applicability  of 
the  method  to  more  complicated  problems.  This  will  be 
explained  later. 

The  Rayleigh-Ritz  method  requires  that  the  number  of  the 
coordinate  functions,  N,  be  determined.  This  is  accomplished 
by  applying  a convergence  criterion  to  the  sequence  of 
approximating  functions  Vj(x),  V2(x),  ...,  v^(x),  where  v^(x) 
(n=l , 2 , . . . , N)  is  a linear  combination  of  coordinate  functions 
as  indicated  in  equation  (2.17).  It  will  be  shown  in  the 
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following  chapters  that  whenever  the  difference  between  the 
Nth  and  (N+l)st  approximation  is  sufficiently  small,  the 
solution  may  be  considered  convergent. 

It  should  be  clearly  understood  that,  in  contrast  to  the 
coefficients  of  a Fourier  expansion,  the  parameters  a^  are 
different  for  each  different  approximation  and  must  be 
re-determined  each  time  a new  approximation  is  chosen [5]. 
However,  a significant  characteristic  that  this  method 
displays  should  also  be  noted.  If 


is  the  Nth  approximation  to  the  true  function  v(x)  and 


is  the  (N+l)st  approximation,  then  the  latter  will  be  better, 
or  at  least  no  worse,  than  the  former.  By  taking  a large 
enough  number  of  the  coordinate  functions,  a convergent  result 
can  be  obtained. 

The  seguence  of  approximations  of  the  Rayleigh-Ritz 
method  will  converge  in  a mean  square  sense  throughout  the 
region  of  interest  (Xo,x^)[5],  i.e.. 


(2.20) 


N+1  * 


(2.21) 


[v(x)  - S a^$Jx)  dx  < e 

n=i  " " 


(2.22) 


where  e is  a positive  real  number,  and 
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■»xl 

[v(x)  - dx  = 0 (2.23) 

xO 

provided  that  the  set  of  coordinate  functions  is  complete. 
An  example  of  a set  of  complete  functions  is  the  complete 
Fourier  series,  including  both  sine  and  cosine  functions. 
This  series  will  be  used  as  the  coordinate  functions  in  the 
following  investigation.  However,  other  functions  can  also 
be  chosen  to  serve  as  coordinate  functions. 

The  Rayleigh-Ritz  method  can  also  be  applied  to  problems 
with  two  or  more  independent  variables,  such  as  problems  of 
plates  and  shells.  For  the  case  of  two  independent  variables, 
the  approximation  can  be  written  as 

v(x,y)  « v^(x,y)  = a^#Jx,y)  (2.24) 

or 

Nm  Nn 

v(x,y)  « v,(x,y)  = a,,f  Jx)  0,  (y)  (2.25) 

where  a^  or  ajj,,,  are  parameters  to  be  determined  such  that  H(aj,) 
or  n(a^^)  is  stationary. 

Laqranaian  Multiplier  Method 

Consider  the  problem  of  extremizing  a function  of  N 


variables  a^,  a^,  ...,  a^,: 

n — II(a^,a2,  ... 

(2.26) 

subject  to  the  constraints: 

i— 1 ,2,...,!^ 

(2.27) 

lim 

n-^oo 
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where  1^,  < N must  hold.  Since  II  is  made  an  extremum, 

TT  N sn 

dll  = E,  — da_  = 0 (2.28) 

n=i  3a 

n 

Because  of  the  constraint  equations  (2.27),  not  all  a^  are 
independent.  Only  (N-I^,)  of  the  variables  a^  are  independent. 
Therefore,  the  differentials  da^  are  not  arbitrary  and  are 
subject  to  the  constraints: 

N 3Gj 

dG^  = S ^da^  = 0 i=l,2,...,I^  (2.29) 

n=i  da^ 

Multiplying  equations  (2.29)  with  the  Lagrangian  multipliers 
(i=l,2,  . . . ,lj,)  and  adding  them  to  equation  (2.28)  yields 

dll  + S A.dG, 

N 3n  ic  3G. 

= S(—  + S A^-^)da,  (2.30) 

n=i  da^  1=1  da^ 

Since  are  arbitrary,  it  is  always  possible  to  choose  them 
in  such  a way  that  of  the  N parenthesized  terms  will 
vanish,  and  only  (N-I^,)  sets  of  terms  are  left.  Furthermore, 
since  the  remaining  (N-I^,)  differentials,  da^,  can  be 

considered  independent,  the  terms  inside  the  remaining 
parentheses  must  be  zero.  Therefore,  the  stationary  point  can 
be  found  directly  from 

311  Ic  3G. 

— — + ,S  A^-^^  = 0 n=l,2,...,N  (2.31) 

3ap  1=1  3a, 


together  with  equations  (2.27).  Equations  (2.27)  and  (2.31) 
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consist  of  (N+Ij.)  equations  with  (N+I^,)  unknowns 

• * • /^Ic)  * 

The  Lagrangian  function,  L,  is  formed  as  follows: 

L ( aj , a2 , . . . , aN , Ai , ^2 / • • • / ^ic)  ~ 

Ic 

n ( 3-2  / / • • • / 3^ ) A^Gj(3j,  32  / •••/3j^)  (2*32) 


The  stationary  point  of  the  function  L can  be  obtained  by 
taking  the  partial  derivatives  of  L with  respect  to  the  (N+I^,) 
variables  aj,  a2,  a^,  Aj,  •••/  ^ic  equating  them  to 

zero: 


3L  / 3a^  = 0 n=l,2,  . . . ,N 


(2.33) 


dL  / dX^  = 0 i=l,2,  . . . ,1^ 

Performing  the  differentiations  of  the  above  equations  also 
results  in  equations  (2.27)  and  (2.31).  This  is  the 
Lagrangian  rule[4][32],  and  it  provides  a concise  way  to 
simultaneously  obtain  the  optimization  equation  (2.31)  and  the 
constraint  equation  (2.27). 

By  using  the  Lagrangian  rule  in  the  Rayleigh-Ritz  method, 
the  satisfaction  of  the  boundary  conditions  required  by  the 
Rayleigh-Ritz  method  is  no  longer  necessary.  The 
applicability  of  the  Rayleigh-Ritz  method  can,  thus,  be 
extended  to  more  complicated  problems  with  a wide  variety  of 
boundary  conditions.  This  will  be  explained  in  the  next 


section. 


Ravleiah-Ritz  Method  Using 
Laaranaian  Multipliers 
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Conventionally,  the  coordinate  functions  ^p(x)  of  the 
Rayleigh-Ritz  method  must  be  so  chosen  that  the  boundary 
conditions  are  all  satisfied.  The  Lagrangian  multiplier 
method  provides  a measure  to  overcome  this  limitation.  By 
applying  the  Lagrangian  multiplier  method,  it  is  not  necessary 
to  limit  the  choice  of  coordinate  functions  by  reguiring  that 
the  boundary  conditions  be  satisfied.  The  unsatisfied 
boundary  conditions  may  be  treated  as  constraint  equations 
(2.27) [18] [19] [20] . 

Since  the  configuration  of  a deformed  body  can  be 
approximated  by  a linear  combination  of  the  coordinate 
functions  with  the  unknown  parameters  a,^,  the  boundary 
conditions  can  also  be  expressed  as  a combination  of  these 
coordinate  functions  and,  thus,  as  a function  of  the  unknown 
parameters.  If  the  unsatisfied  boundary  conditions  are 
regarded  as  the  constraint  equations,  then  the  problem  of 
extremizing  the  potential  function  II  becomes  a constrained 
optimization  problem.  The  number  of  constraints  is  dependent 
upon  how  many  boundary  conditions  are  not  satisfied.  By 
applying  the  Lagrangian  method,  the  same  equation  as  equation 
(2.32)  can  be  obtained  with  II  denoting  the  total  potential 
energy  of  a given  elastic  system  and  denoting  the 
unsatisfied  boundary  conditions.  The  unknown  parameters  a^. 
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as  well  as  the  Lagrangian  multipliers  can  then  be  easily 
obtained. 

It  is  also  possible  to  use  the  same  approach  to  transfer 
any  geometrical  limitation  on  the  deformation  of  the  body  to 
a constraint  eguation,  such  as  in  the  case  of  a multi-span 
beam  where  the  interior  supports  restrict  the  deformation  of 
the  beam.  The  restrictions  caused  by  the  interior  supports 
may  be  regarded  as  constraints,  and  transferred  into  eguations 
that  can  be  expressed  in  terms  of  the  parameters  a,^.  These 
constraints  can  be  included  in  the  process  of  extremizing  the 
potential  energy,  and  the  solution  so  obtained  will  satisfy 
all  boundary  conditions  together  with  those  geometrical 
constraints  imposed  on  the  interior  of  the  deformable  body. 

The  use  of  Lagrangian  multipliers  to  modify  the  Rayleigh- 
Ritz  method  is  the  approach  that  is  proposed  in  this  study  to 
handle  the  analysis  and  optimum  design  of  elastic  structures. 
In  order  to  maintain  focus  on  the  essence  of  this  new 
approach,  only  simple  structural  elements  such  as  beams, 
plates,  and  shells  with  linearly  elastic  deformations  will  be 
considered.  The  principle  of  superposition  will  therefore 
hold  for  all  the  cases  studied.  The  illustrated  examples, 
although  simple,  exhibit  all  the  basic  features  of  more 
complex  problems.  The  next  chapter  will  be  devoted 
exclusively  to  beam  problems  where  the  application  of  the 
method  can  be  illustrated  with  minimum  complexity. 


CHAPTER  3 

APPLICATIONS  TO  BEAM  PROBLEMS 


Introduction 

The  governing  equations  of  an  elastostatic  problem  can 
be  grouped  into  three  sets [33]: 

1.  Equilibrium  equations  in  terms  of  stress  components 

2.  Compatibility  equations  in  terms  of  strain  components 

3 . Constitutive  equations  relating  the  stresses  and  the 
strains 

Boundary  conditions  corresponding  to  set  1 are  named  force 
boundary  conditions  and  those  corresponding  to  set  2 are 
called  geometry  boundary  conditions.  All  boundary  conditions 
fall  into  these  two  categories. 

For  a beam  subjected  to  a lateral  load  p(x),  the 
equilibrium  equation  is [15] 


d^M(x) 

— 71 — = -P(x) 
dx 


(3.1) 


where  M(x)  denotes  the  bending  moment  of  the  beam  and  x is  the 
coordinate  along  the  beam.  The  compatibility  equation  is [15] 


d^v  (x) 
dx^ 


= -p{x) 


(3.2) 


where  v(x)  is  the  deflection  function  of  the  beam  and  fi (x)  is 
the  curvature.  The  constitutive  equation  requires  that [15] 
M(x)  = -EI/)(x)  (3.3) 
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By  combining  these  equations,  the  governing  differential 
equation  of  a beam  subjected  to  a lateral  load  is  obtained. 
This  was  shown  in  equation  (2.16).  Since  the  differential 
equation  is  of  the  fourth  order,  four  boundary  conditions  are 
necessary  to  complete  the  solution.  Usually,  the  boundary 
conditions  at  each  end  of  the  beam  are  a combination  of  any 
two  of  the  following: 

i.  v(x) 

ii.  dv(x)/dx 

iii.  M(x) 

iv.  V(x)  = dM(x)/dx 

where  V(x)  denotes  the  shear  force.  Equations  (i)  and  (ii) 
are  geometry  boundary  conditions  which  relate  to  generalized 
displacements  of  the  deflection  of  the  beam,  and  the  remaining 
two  are  force  boundary  conditions  which  relate,  as  the  name 
suggests,  to  generalized  forces.  Using  equation  (3.2), 
equations  (iii)  and  (iv)  can  be  rewritten  as 


It  is  therefore  possible  to  express  any  boundary  conditions 
in  terms  of  the  deflection  function  v(x) . Since  only  linearly 
ealstic  deformations  are  considered  in  this  study,  the  above 
expressions  will  be  applied  only  to  beams  deflecting  in  the 
linearly  elastic  range. 


M(x) 


(3.4.1) 


and 


V(x) 


(3.4.2) 
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In  analyzing  beam  problems  with  the  Rayleigh-Ritz  method, 
the  deflection  of  the  beam  is  assumed  to  be  a combination  of 
some  suitably  chosen  coordinate  functions.  For  all  the 
examples  discussed  below,  Fourier  series  are  chosen  as  the 
coordinate  functions  because  they  satisfy  a wide  range  of 
boundary  conditions.  If  the  coordinate  system  is  taken  the 
same  as  that  shown  in  figure  3.1,  then  the  deflection  function 
can  be  assumed  as 


® nrrx 

v(x)  =Sja,sin(-— ) 


(3.5) 


where  I is  the  length  of  the  beam.  If  the  bending  rigidity 
El  is  constant,  the  strain  energy  of  the  beam  is [15] 


U = 


El  d^v 


EItt^ 


( ? ) dx  = o _ 

2 dx^  4£^  n=l  " 


2 a„  n^ 


(3.6) 


The  work  done  by  the  external  load  p(x)  acting  on  the  beam 
takes  the  following  form: 


W = 


p(x)v(x)dx  = S a^ 

n=l  " 


0 


nTTX 

p(x)sin(— ^)dx  (3.7) 


Accordingly,  the  total  potential  energy  of  the  beam  is 
n = U - W 


EItt  •>  0 a ® 

2 a„  n - S a„ 


4 £3  n=l  " 


n=l 


0 


nTTX 

p(x)sin( — ^)dx  (3.8) 


Simplv-Supported  Beams 


For  the  simply-supported  beam  shown  in  figure  3.1,  the 
boundary  conditions  for  both  ends  are  v=0  and  M=0  (or  v”=0. 
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if  divided  by  the  bending  rigidity  El) . Since  sine  functions 
already  satisfy  these  boundary  conditions  at  both  x=0  and  x=£ , 
no  constraint  is  necessary  for  this  case. 


p(x) 


v(x) 


X 


V 


Figure  3.1  A simply-supported  beam 


If  p(x)  is  a uniform  load,  i.e.,  p(x)=Pq,  then 


W = S a 


Po^ 

n=l  nTT 


(1-COS  (nTT)  ) 


and 


(3.9) 


n = 


■H!L.  ^ _ 

4^3  n=l  n 


Pf)£  ® 1 - cos  (nTT) 

— — S a„ 

n=l  " 


7T 


n 


(3.10) 


Taking  partial  derivatives  of  equation  (3.10)  with  respect  to 
the  parameters  a^  and  setting  them  equal  to  zero  gives 
3n  ElTr'^n'*  Po^  ( i ” cos  (nTT)  ) 


3a„ 


2Z 


3 “n 


a„  - 


nTT 


= 0 (3.11) 


Thus  the  parameters  a^  can  be  obtained  as 


= 


2Pq£  ( 1 - cos(nTr)  ) 

ElnV® 


(3.12) 


It  should  be  noted  that  in  this  case  these  parameters  are  also 
the  coefficients  of  the  Fourier  expansion  of  the  theoretical 
deflection  of  a simply-supported  beam  subjected  to  a uniformly 
distributed  load[15] [34] . The  Nth  approximation 
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N nTTX 

v^(x)  =S^a,sin(  — ) (3.13) 

is  simply  a truncated  Fourier  series,  and  it  gets  closer  to 
v(x)  as  N becomes  larger.  Although  it  happens  that  the 
parameters  a^^  are  equal  to  the  Fourier  coefficients  for  this 
case,  in  general  they  are  different.  When  the  boundary 
conditions  are  not  identically  satisfied,  the  parameters  a,^ 
deviate  from  the  Fourier  coefficients.  This  can  be 
illustrated  by  the  following  case. 

Clamped-Supported  Beams 

Consider  a clamped-supported  beam  acted  upon  by  a 
uniformly  distributed  load  p(x)=Pq  as  shown  in  figure  3.2. 


V 

Figure  3.2  A clamped-supported  beam 

By  choosing  the  same  coordinate  functions,  the  formulation  of 
the  total  potential  energy  is  the  same  as  that  of  the  simply- 
supported  beam  case.  The  boundary  conditions,  however,  are 
not  all  satisfied  by  the  coordinate  functions.  The 
unsatisfied  boundary  condition  is  v'=0  at  x=0,  which  can  be 
written  as 
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Gi(ap)  = 2^  a^n  = 0 


(3.14) 


Applying  the  proposed  modification  to  the  Rayleigh-Ritz 
method,  the  Lagrangian  function  takes  the  form 
L = n + lGi(aJ 


, . Ppf  <0  1 - cos(nTr)  « 

" ^ — + \ S an  (3.15) 

n=l  ' 


EItt^  <0 

a_‘n‘*  - — S a„ 

n=l  " ^ n=l 


n 


By  taking  partial  derivatives  of  the  Lagrangian  with  respect 
to  a„  as  before,  it  is  found  that 


2Z 


= 


Pnf  ( 1 - cos(nTr)  ) 

4 4 ( An  ) (3.16) 

EItt  n n7T 

Substituting  a^  into  equation  (3.14),  the  undetermined 

multiplier  A can  be  found  by  the  following  equation: 


As  — o = S 

n=l  ^=1 


® Po^  ( 1 - cos(nTr)  ) 


n n 


(3.17) 


The  calculation  of  the  undetermined  multiplier  A is  generally 
not  possible  if  the  upper  limit  of  the  summation  is  infinity. 
To  make  it  feasible  to  calculate  A,  the  first  N terms  of  the 
series  (3.5)  are  taken  as  an  approximation  (the  Nth 
approximation  as  described  in  Chapter  2)  . In  this  case, 
equation  (3.17)  becomes 


N 1 N P(,£(  1 - COS(n7T)  ) 
A 2j  o — ^ 


n=l  n" 


n=l 


n TT 


(3.18) 


The  undetermined  multiplier  A can  now  be  obtained  easily.  By 
substituting  A in  equation  (3.16),  the  parameters  a^  can  be 
calculated  and  thus  the  deflection  of  the  beam  can  be  found. 
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Table  3.1  gives  a comparison  of  the  Fourier  coefficients  of 
the  theoretical  deflection  function  with  the  parameters  a^ 
from  the  modified  Rayleigh— Ritz  method.  As  can  be  seen  from 
the  table,  the  parameters  of  the  new  method  get  closer  to  the 
Fourier  coefficient  as  N gets  larger. 


Table  3.1  Comparison  of  the  Fourier  series  with 
three  approximations 


n 

F.  S. 

5th  App. 

10th  App. 

20th  App. 

1 

5.00817 

4.01582 

4.51443 

4.76198 

2 

-1.00786 

-1.13190 

-1.06958 

-1.03863 

3 

-0.24484 

-0.28159 

-0.26312 

-0.25395 

4 

-0.12598 

-0.14149 

-0.13370 

-0.12983 

5 

-0.06032 

-0.06826 

-0.06427 

-0.06229 

Figure  3.3  shows  the  deflection  curves  of  the  three 
different  approximations  indicated  in  table  3.1,  together  with 
the  theoretical  deflection  curve [15].  The  deflections  shown 
in  the  figure  have  been  normalized  by  the  expression  Pq^Vei. 
It  can  be  seen  from  figure  3.3  that  as  N increases,  the 
approximations  approach  the  theoretical  result.  The  error  at 
x=£/2  is  about  4.6%  for  the  20th  approximation. 

Clamped-Clamped  Beams 

One  of  the  advantages  of  the  proposed  method  is  that  it 
is  easy  to  modify  the  boundary  conditions  without  any  major 
change  in  the  formulation.  Taking  the  previous  case  as  an 
example,  it  is  easy  to  obtain  the  deflection  curve  of  a 
clamped-clamped  beam  (figure  3.4)  from  the  previous 
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Figure  3.3  Deflection  curves  of  a clamped-supported  beam 
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formulation  by  simply  adding  one  more  constraint  to  the 
supported  end  of  the  clamped-supported  beam.  This  forces  the 
slope  of  the  deflection  curve  to  be  zero  at  the  supported  end. 
If  the  Nth  approximation  is  taken,  then  this  constraint  is 

N 


G2(a^)  = a^n  cos(nTr)  = 0 


(3.19) 


Since  the  approximate  function  of  (3.13)  has  already 
satisfied  the  requirement  of  zero  deflection  at  both  ends  of 
the  beam,  equation  (3.19)  and  equation  (3.14)  constitute  the 
only  constraints  needed  for  the  present  case. 


Figure  3.4  A clamped-clamped  beam 

The  Lagrangian  function  takes  the  form 

L = n + (3.20) 

where  denotes  the  number  of  constraints  to  be  imposed  on 
the  beam,  and  is  two  in  this  case  (the  approximate  deflection 
function  v,g(x)  of  equation  (3.13)  has  already  satisfied  the 
boundary  conditions  of  zero  deflection  at  both  ends  of  the 
beam) . Following  the  same  procedure  used  in  the  previous 
section,  it  is  found  that 
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Po^  ( 1 “ cos(n7T)  ) 
ElTr'^n'^  ^ n7T 


AjH  - >211  cos(nTr)] 
(3.21) 


By  substituting  into  Gj  and  G2,  two  algebraic  equations  with 
two  unknowns  Aj  and  >2  are  obtained.  By  solving  these 
equations  and  substituting  the  result  in  equation  (3.21),  the 
displacement  parameters  a^  can  be  obtained.  The  results  for 
the  10th,  the  30th,  and  the  50th  approximations  are  plotted 
in  figure  3.5.  The  deflections  shown  in  the  figure  have  been 
normalized  by  the  expression  Pq^Vei.  The  theoretical 
def lection [ 15]  is  also  plotted  in  the  same  figure  for 
comparison.  The  error  of  the  30th  approximation  at  the  center 
of  the  beam  is  about  5.5%  and  that  for  the  50th  approximation 
is  3.3%. 


Cantilever  Beams 


Under  certain  conditions,  for  example  a cantilever  beam, 
it  is  necessary  to  use  more  complicated  series,  such  as  ones 
including  both  sine  and  cosine  functions,  as  the  coordinate 
functions.  The  nature  of  the  sine  functions  alone  makes  it 
impossible  to  provide  a deflection  at  the  free  end  and, 
similarly,  it  is  not  possible  to  obtain  a slope  other  than 
zero  at  the  free  end  from  the  cosine  functions  alone. 
Consider  the  cantilever  beam  shown  in  figure  3.6.  A uniform 
load  Pg  is  distributed  over  the  left  half  of  the  beam. 

Assume  that  the  defection  curve  is  a linear  combination 
of  both  sine  and  cosine  functions  as  follows: 
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Figure  3.5  Deflection  curves  of  a clamped-clamped  beam 
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Vn(x) 


N 

3q  + 2^  a^cos  (n7TX/£) 


N 

+ S b sin(n7rx/£) 

n=l  ' 


(3.22) 


The  strain  energy  stored  in  the  beam  for  this  assumed 
deflection  is 


U 


EItt'* 


N 

2 

n=l 


(an'+bn')n' 


+ 


N N 4n^m^ 

n?i  ■.  2_^2.  ( COS  ( n7T ) COS  (mTT ) " 1 ) ] (3.23) 

nff  n 7t  {ri  in  ; 


Figure  3.6  A cantilever  beam 


The  work  done  by  the  external  load  is 


W 


Po^Qg 

2 


N 

+ 2 
n=l 


Po^n^ 

nn 


nn 

sin(  — ) 
2 


N 

+ 2 
n=l 


Ppbn^ 

nn 


nn 

(l-cos(  — ) ) 

(3.24) 


The  boundary  conditions  are  v=0  and  v*=0  at  x=0,  and  v”=0  and 
v"'=0  at  x=£ , none  of  which  are  satisfied  by  the  approximation 
(3.22).  The  constraints  are  thus 


Gi  = 


N 

+ 2 
n=l 


= 0 


(3.25.1) 


N 

Go  = 2 b„n  = 

n=l 


0 


(3.25.2) 
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N 2 

G3  = a^n  COS  (n7T)  =0  (3.25.3) 

N 3 

G4  = b^n  COS  (n7T)  =0  (3.25.4) 

The  Lagrangian  function  takes  the  same  form  as  equation  (3.20) 
with  Ij,=4 . By  taking  partial  derivatives  of  the  Lagrangian 
with  respect  to  parameters  a^,  a^,  and  (n=l,  2 , . . . ,N)  , a set 
of  (2N+1)  linear  algebraic  equations  with  ag,  a^^,  b^ 
(n=l,  2 , . . . ,N)  , and  (i=l,  2 , . . . , 1^,)  as  unknowns  is  obtained. 
These  equations  take  the  following  form: 

EItt'*  4 N 4n^m^ 

^ [ 2a„n  + S b^ 5 — 5- (cos  (n7T)  cos  (mTT) -1)  ] 

n "’Tr(nW)' 

^ ^ 2 Po^  . 

+ A,  + A3ncos(n7T)  = — ^sin( ) 

nir  2 

n=l,2,...,N  (3.26.1) 

EItt^  4 N 4m^n^ 

— ^ [ 2b^n  + E /a^— — 2 — ^ (cos  (nTr)  cos  (m7r) -1)  ] 

4f  ni^iv  7r(m"-n‘') 

, ^ 3 Pgf  rnr 

+ Agn  + A4ncos(n7T)  = — 1-cos  ( )] 

n7T  2 

n=l,2,...,N  (3.26.2) 

^1  " Po^/2  (3.26.3) 

The  parameters  ag,  a^,  and  b^  (n=l,2,  . . . ,N)  , and  the  Lagrangian 
multipliers,  A^  (i=l, 2 , . . . , 1^,)  , can  be  obtained  by  solving  the 
above  equations,  together  with  the  constraint  equations 
(3.25).  Figure  3.7  shows  the  comparison  of  the  fifth  and  the 
tenth  approximations  with  the  theoretical  solution [ 15 ] . The 
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deflections  shown  in  the  figure  have  been  normalized  by  the 
expression  Pq^Vei. 

As  can  be  seen  from  figure  3.7,  the  10th  approximation 
almost  coincides  with  the  theoretical  curve.  A very  good 
approximation  can  be  obtained  by  using  only  ten  cosine  and 
sine  functions.  Again,  it  should  be  noted  that  using  only 
sine  functions  or  only  cosine  functions  would  not  result  in 
an  acceptable  solution. 

It  is  also  possible  to  obtain  a good  approximation  by 
replacing  the  cosine  functions  with  a linear  term,  i.e., 
Vo(x)=aoX.  To  do  so,  Vq(x)  should  be  added  to  equation  (3.13), 
and  then  the  same  procedure  followed.  Since  the  second 
derivative  of  Vg(x)  vanishes,  the  formulation  for  strain 
energy  is  not  affected.  Only  the  formulation  for  work  would 
need  to  be  modified. 

Cantilever  Beams  with  the  Free  End  Flexibly  Supported 

As  shown  in  the  previous  case,  the  method  presented  in 
this  investigation  makes  it  very  easy  to  change  the  boundary 
conditions  of  a beam  problem  without  any  major  modification 
of  the  formulation.  To  further  extend  this  method,  let  us 
consider  a cantilever  beam  with  the  free  end  supported  by  a 
spring  with  stiffness  K as  shown  in  figure  3.8.  Although  in 
the  classical  method,  this  type  of  problem  is  treated 
differently  from  the  previous  case,  the  proposed  method 
produces  a formulation  that  is  quite  similar. 
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Figure  3 . 7 Deflection  curves  of  a cantilever  beam 
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V 


Figure  3.8  A cantilever  beam  with  the  free  end  flexibly 
supported 

The  spring  stores  energy  when  it  is  deflected.  The 
amount  of  potential  energy  stored  in  the  spring  should  be 
added  to  the  total  potential  energy.  Assuming  the  spring  is 
not  preloaded  and  the  deflection  of  the  spring  is  rj,  the 
potential  energy  of  the  spring  can  be  written  as 

= Kr}^/2  (3.27) 
Compatibility  between  the  compression  of  the  spring  and  the 
deflection  of  the  beam  requires  that 

N 

Gg  = ( ao  a^cos(nTr)  ) - 77  = 0 (3.28) 

This  equation  is  also  an  equality  constraint  equation. 

Following  the  same  procedure  as  in  the  previous  section, 
similar  equations  can  be  obtained,  along  with  one  additional 
equation  obtained  by  taking  the  partial  derivative  of  L with 
respect  to  rj,  which  is 

Kt]  - A5  = 0 (3.29) 
where  Ag  is  the  undetermined  multiplier  associated  with  the 
constraint  equation  Gg.  It  can  be  seen  from  equation  (3.29) 
that  the  Lagrangian  multiplier  associated  with  the  constraint 
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of  the  deflection  of  the  beam  and  the  compression  of  the 
spring  is  the  reaction  force  between  them.  In  this  situation, 
the  Lagrangian  multiplier  associated  with  the  compression  of 
the  spring  has  the  dimension  of  force.  The  product  of  the 
multiplier  and  the  deflection  constraint,  which  has  the 
dimension  of  deflection,  has  the  dimension  of  work.  This 
makes  the  dimension  of  the  Lagrangian  function  consistent. 
It  may  be  deduced  from  these  observations  that  multipliers 
associated  with  generalized  displacement  constraints 
(deflections  or  slopes  of  the  deflection  curve)  have  the 
dimensions  of  generalized  forces.  Furthermore,  it  can  also 
be  deduced  that  the  multipliers  are  the  generalized  forces 
required  to  insure  compliance  with  the  desired  boundary 
conditions.  In  other  words,  forces  are  superimposed  on  the 
elastic  body  through  the  use  of  the  undetermined  multipliers 
to  modify  the  boundary  conditions  as  required. 

Solving  equations  (3.25),  (3.26),  (3.28),  and  (3.29) 
simultaneously  can  give  the  deflection  of  the  beam.  Taking 
K=18EI/£^,  figure  3.9  shows  the  result  obtained  from  the 
classical  method  and  the  result  of  the  fifth  approximation. 
The  deflections  shown  in  the  figure  have  been  normalized  by 
the  expression  It  can  be  seen  that  the  error  is 
negligible  in  the  plot.  The  compression  of  the  spring 
obtained  from  the  proposed  method  is  2 . 60682x10’^  (Pq^Vei)  and 
that  obtained  from  the  classical  method  is  (1/384)  (Pq^VeI)  ; 
the  error  is  approximately  1%.  The  reaction  force  obtained 
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Figure  3.9  Deflection  curves  of  a cantilever  beam  with  the 
free  end  flexibly  supported 
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from  the  proposed  method  is  4 . 69228xlO’^Pg£ , which  is  about 
0.1%  more  than  the  value  calculated  from  the  classical  method. 

Beams  Supported  on  Elastic  Foundations 

By  extending  the  previous  case,  a beam  placed  on  an 
elastic  foundation  can  be  treated  in  a similar  way.  The 
elastic  foundation  can  be  regarded  as  a combination  of  a 
number  of  springs.  Consider  the  following  situation  (figure 
3.10),  a simply-supported  beam,  resting  on  an  elastic 
foundation  with  modulus  K^,  is  bent  by  a moment  Mq,  acting  at 
one  end  of  the  beam. 


Figure  3.10  A beam  supported  on  an  elastic  foundation 

By  using  the  approximation  series  (3.13) , the  formulation 
of  the  strain  energy  of  the  beam  is  the  same  as  equation  (3.6) 
except  that  the  upper  limit  of  the  summation  is  replaced  by 
N.  The  energy  stored  in  the  foundation  is 


Mo 


7^/  / / y / y 


i 


X 


V 


Uf 


2 

J 0 


(3.30) 


The  work  done  by  Mg  is 


W 


(3.31) 
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The  Lagrangian  function  takes  the  following  form: 


L = 


EItt^  n 


2_4 


n=l 


S a/n^  + 


2 a 

n=l  " 


M(,7T  N 

— ^ 2 an 

P n=l  " 


(3.32) 


By  taking  partial  derivatives  of  L with  respect  to  a^  and 
solving  for  a^,  it  is  found  that 


MgTrn 


/ ( 


ElTr'^n'^ 

21^ 


+ 


) 


(3.33) 


Figure  3.11  shows  the  comparison  between  the  theoretical 
solution  [35]  and  the  results  of  the  third  and  the  fifth 
approximations,  in  which  is  taken  as  64EI.  The  deflections 
shown  in  the  figure  have  been  normalized  by  the  expression 
Mq^Vei.  The  figure  shows  that  only  five  sine  functions  give 
an  excellent  fit  to  the  theoretic  curve. 


Multi-Span  Beams 

It  has  been  mentioned  previously  that  the  proposed  method 
can  be  applied  to  solve  the  problem  of  multi-span  beams.  In 
this  case,  a two-span  beam  (figure  3.12)  is  taken  as  an 
example  to  demonstrate  the  applicability  of  the  method.  By 
taking  the  same  approximation  series  as  eguation  (3.22),  the 
formulation  for  strain  energy  of  the  beam  takes  the  same  form 
as  eguation  (3.23).  The  work  done  by  the  concentrated  load 
P is 

W = P v^(Xp) 

= P [ ag  + 2^  a^cos(n7TXp/£)  + b^sin  (n7rXp/£ ) ] 


(3.34) 
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Figure  3.11  Deflection  curves  of  a beam  supported  on  an 
elastic  foundation 
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The  constraint  equations  due  to  the  unsatisfied  boundary 
conditions  are  given  as 

Gi  = ao  + a,  = 0 (3.35.1) 

G2=Jia,n^  = 0 (3.35.2) 

N 

^ ^0  + (”’'■)  = 0 (3.35.3) 

N 2 

G4  = a^n  cos(n7T)  = 0 (3.35.4) 


P 


Xp 


Figure  3.12  A multi-span  beam 


The  interior  support  will  prevent  the  beam  from  bending 
under  the  effect  of  the  force  P and,  thus,  will  constrain  the 
deflection  of  the  beam  at  x=x^.  This  gives  an  additional 
constraint  equation  which  must  be  satisfied  as  well  as  the 
boundary  conditions  in  the  process  of  extremizing  the 
potential  energy.  This  constraint  can  be  put  in  the  following 


form: 
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^5  = ^0  + a^cosCnTTXg/^)  + b^sin(n7rx3/£)  = ° 

(3.35.5) 

Using  the  same  procedure  of  solving  for  the  parameters  a^,  a^, 
and  (n=l,  2 , . . . ,N)  as  in  the  previous  cases,  the  result 
shown  in  figure  3.13  is  obtained.  Taking  N=5  gives  an 
approximation  having  almost  no  difference  from  the  theoretical 
result  obtained  with  classical  techniques.  The  classical 
method  usually  takes  two  steps;  first  solving  the  deflection 
due  to  the  external  force  P and  determining  the  reaction  of 
the  interior  support,  and  then  superimposing  the  deflections 
due  to  both  force  P and  the  reaction.  If  there  is  more  than 
one  interior  support,  the  classical  method  becomes  tedious. 
However,  for  the  proposed  method,  the  only  difference  is  the 
number  of  constraints.  This  does  not  change  the  basic 
formulation  and  does  not  introduce  any  additional  complexity 
to  the  solution  process. 
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Figure  3.13  Deflection  curves  of  a multi-span  beam 


CHAPTER  4 

APPLICATIONS  TO  PLATE  PROBLEMS 
Introduction 

In  this  investigation,  a plate  is  defined  as  an  elastic 
body  bounded  by  two  flat  surfaces.  The  distance  between  the 
surfaces  is  the  thickness  of  the  plate.  The  plane  parallel 
to  the  surfaces  and  bisecting  the  thickness  is  called  the 
middle  plane.  The  bending  property  of  the  plate  depends  upon 
its  dimensions,  especially  the  thickness.  Based  on  the  amount 
of  deflection  and  the  thickness  of  the  plate,  plate  problems 
can  be  grouped  into  three  categories [ 3 ] : 

1.  thin  plate  with  small  deflection 

2 . thin  plate  with  large  deflection 

3.  thick  plate 

Usually,  thick  plate  problems  are  treated  as  three-dimensional 
problems.  The  problem  of  a thin  plate  with  large  deflections 
is  a nonlinear  elasticity  problem.  Both  of  these  categories 
are  not  considered  in  this  investigation.  The  first  category, 
which  is  the  only  one  to  be  dealt  with  in  this  study,  is 
considered  linearly  elastic  in  classical  plate  theory. 

The  development  of  the  theory  of  bending  of  a thin  plate 
is  based  on  the  following  assumptions  attributed  to 
Kirchhoff [3] [16] [36] : 
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1.  The  linear  filament  of  the  plate  which  is  initially 
normal  to  the  middle  plane  remains  straight  and 
normal  to  the  middle  plane  after  bending. 

2.  There  is  no  deformation  in  the  middle  plane,  i.e., 
this  plane  remains  neutral  and  unstretched. 

3 . Normal  stress  in  the  direction  transverse  to  the 
middle  plane  is  small  compared  to  the  other  normal 
stresses,  and  the  normal  strain  in  that  direction  is 
small  enough  to  be  disregarded. 

These  assumptions  are  good  if  the  deflection  is  small  in 

comparison  with  the  thickness.  They  enable  us  to  express  the 

potential  energy  of  the  bent  plate  in  terms  of  the  curvatures 

of  the  bent  plate.  Furthermore,  they  also  make  it  possible 

to  express  all  the  stresses  in  terms  of  the  deflection  of  the 

plate. 

In  the  following  discussion,  the  middle  plane  is  supposed 
to  be  coincident  with  the  xy  plane  and  the  deflection  is  in 
the  z direction.  Thus,  the  deflection  is  a function  of  both 
X and  y.  With  this  in  mind,  the  compatibility  eguation 
requires  that 

/>x  = = -d^v/dx^  (4.1) 

and  similarly, 

Py  = 1/ry  = -d^w/dy^  (4.2) 

where  p denotes  the  curvatures,  r denotes  the  corresponding 
radii  and  v is  the  deflection  function  of  the  plate.  The 
subscript  x denotes  the  x-direction  and  y denotes  the  y- 
direction.  The  equilibrium  of  forces  requires  that [3] 

+ a^My/ay^  - a^M^y/axSy  = -p(x,y)  (4.3) 

where  and  My  denote,  respectively,  the  bending  moment  in  the 
X and  y directions.  The  symbol  M,^y  denotes  the  twisting  moment 
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and  p(x,y)  is  the  transverse  load  acting  upon  a surface  of  the 
plate. 

The  constitutive  equation  requires  that [3] 

= D{p^+vp^)  = -D(3^v/3x^+J/a^v/3y^)  (4.4) 

= D(p^+j//0y)  = -D(a^v/ax^+j/aV/3y^)  (4.5) 

M^y  = D(l-l/)  aV/ax3y  (4.6) 

where  v is  Poisson's  ratio.  The  symbol  D is  the  flexural 
rigidity  of  the  plate  and  is  defined  as 

D = EhV[12(l-J/)  ] (4.7) 


In  equation  (4.7)  , E is  Young's  modulus  and  h is  the  thickness 
of  the  plate. 

By  combining  these  equations,  the  governing  equation  of 
the  plate  is  obtained  as  follows [3]; 


a V 


+ 2 


aV 


+ 


a V 


D 


(4.8) 


ax"  ax^ay"  ay" 

Equation  (4.8)  is  the  Lagrange  equation  which  relates  the 
transverse  deflection  v(x,y)  to  the  applied  load  p(x,y). 

The  boundary  conditions  at  the  edges  of  the  plate  can  be 
categorized  as  follows[3] [16] [36] : 

1.  clamped  edge  along  which  v and  ^v/^T^  vanish 

2.  simply-supported  edge  along  which  v and  vanish 

3 . free  edge  along  which  and  vanish 

In  the  above  statements,  the  derivative  with  respect  to  n 
means  a derivative  with  respect  to  the  normal  to  the  edge,  and 
and  denote,  respectively,  the  moment  and  the  shear  force 
in  the  normal  direction  to  the  edges. 
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Equation  (4.8)  , as  well  as  the  above  boundary  conditions, 
can  also  be  obtained  by  the  energy  principle[3] [16] [36] . The 
process  of  variation  was  worked  out  by  Lord  Rayleigh  in  his 
book  "Theory  of  Sound" [37].  It  has  also  been  shown  by  other 
researchers [ 6 ][ 17 ] that  the  solution  to  equation  (4.8) 
minimizes  the  potential  energy  of  the  plate.  The  potential 
energy  takes  the  following  form: 

n = U - W (4.9) 
where  U is  the  strain  energy  of  the  plate,  and  is  obtained  by 
performing  the  integration  of  the  strain  energy  density  of  the 
plate  over  the  whole  plate.  The  strain  energy  U for  a plate 
is  formulated  as [3] 


2(l-u) [ 


3^v  5^v 


dxdy 


)^]  }dA 


(4.10) 


The  symbol  A denotes  the  surface  area  of  the  plate.  If  a load 


P(x,y)  acts  upon  one  surface  of  the  plate,  the  work  done  by  the 
load  is 


W 


P(x,y)v(x,y)  dA 
J A 


(4.11) 


Plates  with  All  Edges  Simply  Supported 


The  problem  of  a plate  with  all  edges  simply  supported 
has  been  cited  by  many  authors,  for  example  Richards [1], 
Timoshenko[3 ] , Saada[16],  and  Love [36],  to  demonstrate  the  usage 
of  the  Rayleigh-Ritz  method  in  plate  problems.  Taking  a rectangular 
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plate  with  dimensions  and  £y  (shown  in  figure  4.1)  as  an  example, 
the  coordinate  functions  are  chosen  to  be  sine  functions  for 
both  the  X and  y directions,  i.e., 

Nm  Nn 

v(x,y)  « S^a^^sin(m7rx/f^)sin(n7ry/£y)  (4.12) 

To  simplify  the  problem,  and  will  be  taken  the  same  throughout 
this  study.  With  the  deflection  surface  approximated  by  the 
series  of  eguation  (4.12),  the  strain  energy  is  expressed  as 


U = 


mTT 


8 m=l  n=l 


riTT 


(4.13) 


which  is  obtained  by  performing  the  integration  of  eguation  (4.10) 


s.s. 


1 

1 

s.s. 

p<x, y) 

1 

S.S, 

1 

s.s. 


Figure  4.1  A plate  with  all  edges  simply  supported 


Similarly,  by  performing  the  integration  of  eguation  (4.11)  , 
the  formula  for  the  work  done  by  the  applied  load  p(x,y)  can 
be  obtained.  This  formula  is  expressed  as  a function  of  the 
parameters  a^.  If  the  external  load  p(x,y)  is  a laniformly  distributed 
load,  i.e.  , p(x,y)=pQ,  then  the  work  done  by  Pq  can  be  formulated 
as 
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W = 


Pn-^y^wNm  Nn 

S ( 

m=l  n=l 


7T 


1-COS  (itlTT)  1-COS  (HTT) 

)( -)  (4.14) 


m 


n 


Since  the  second  derivative  of  a sine  coordinate  function 
vanishes  at  the  edges  of  the  rectangular  plate,  i.e.,  d^v/dx^=0 
at  x=0  and  x=£^  and  d\/dy^=0  at  y=0  and  y=fy,  all  boundary  conditions 
for  the  rectangular  simply— supported  plate  are  satisfied  by  the 
approximation  series  (4.12)  and  no  constraint  is  necessary. 
In  this  case,  the  Lagrangian  function  is  the  potential  energy 
of  the  plate,  i.e.. 


L = n = U - W (4.15) 

Using  the  same  approach  as  in  the  beam  problems,  the  partial 
derivatives  of  L with  respect  to  each  a^^  are  taken,  and  set  equal 
to  zero.  The  displacement  parameters  a^^  are  obtained  as 

^mn  ~ ^mn  / ( ^ ^mn ) (4.16) 

where 


and 


D£  £ mTT 

E = ^ r ( 

mn  « L V 


)'  + ( 


nn 


(4.17) 


P(,£^£^  1-cos  (mTT)  1-cos  (nTT) 

Smn  = ( )( —)  (4.18) 

JT  m n 

In  view  of  the  above  equations,  it  can  be  seen  that  the  numerical 
value  of  a^^  diminishes  very  quickly  as  m and  n increase. 

The  result  for  a square  plate  is  shown  in  figures  4.2  and 
4.3.  Figure  4.2  shows  deflection  curves  taken  along  the  section 
y=^y/2.  The  deflections  have  been  normalized  by  the  expression 
Pq^xVd.  It  can  be  seen  from  figure  4.2  that  for  this  case  it 
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Figure  4.2  Deflection  curves  of  a simply-supported  plate 
taken  along  the  section  y=£  /2 
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is  difficult  to  distinguish  the  three-term  approximation  curve 
from  the  nine-term  approximation  curve.  The  difference  between 
the  one-term  approximation  and  the  nine-term  approximation  is 
about  2.4%  at  the  point  x=i^/2  and  y=£y/2  and  that  between  the 
three-term  approximation  and  nine-term  approximation  is  only 
0.17%.  From  this,  the  sequence  of  approximations  can  be  considered 
convergent  when  and  are  greater  than  three.  Since  the  difference 
between  the  nine-term  approximation  and  the  theoretical  deflection[3] 
is  almost  negligible,  it  can  be  seen  that  taking  only  one  term 
of  the  series  (4.12)  , i.e. , N^=N^=1,  will  give  an  accurate  enough 
result  (with  an  error  of  about  2.4%)  for  engineering  applications. 

Figure  4.3  shows  the  deflection  surface  of  the  square  plate 
obtained  by  the  nine-term  approximation.  The  deflection  in  figure 
4.3  has  been  normalized  in  the  same  manner  as  above,  and  magnified 
25  times  with  respect  to  the  dimensions  of  the  plate. 


Figure  4.3  The  deflection  surface  of 
a simply-supported  plate 


Plates  with  One  Edge  Clamped 
and  Three  Edges  Simply  Supported 
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As  illustrated  in  the  previous  chapter,  the  Rayleigh-Ritz 
method  using  undetermined  multipliers  makes  it  very  convenient 
to  convert  one  beam  problem  to  another  with  different  boundary 
conditions.  It  will  also  be  easy  to  change  the  boundary  conditions 
of  a plate  problem  without  any  major  change  in  the  formulation 
by  using  the  modified  Rayleigh-Ritz  method.  Consider  a rectangular 
plate  with  three  edges  simply  supported  and  the  fourth  edge  clamped 
as  shown  in  figure  4.4. 


s.s. 


X 

Figure  4.4  A plate  with  one  edge  clamped  and  three 
edges  simply  supported 

With  the  same  assumed  deflection  function  as  equation 
(4.12) , the  formulations  of  strain  energy  and  work  will  be  the 
same  as  equations  (4.13)  and  (4.14).  The  boundary  conditions, 
however,  are  not  all  satisfied  by  the  approximate  deflection 
function.  What  is  not  satisfied  is  the  slope  at  the  fixed  boundary, 
i.e. , 


3v  / 3x  = 0 


for  x=0 


(4.19) 
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Since  it  is  extremely  difficult,  if  not  impossible,  to  find  a 
continuous  function  to  serve  as  a "Lagrangian  multiplier"  to 
force  the  slope  to  be  zero  along  the  whole  edge  x=0,  an  alternative 
way  of  applying  the  concept  of  Lagrangian  multipliers  must  be 
taken. 

In  this  alternative  approach,  the  boundary  conditions  are 
only  satisfied  using  undetermined  multipliers  at  a finite  number 
of  points  chosen  from  the  edge  x=0,  and  the  segments  in  between 
are  disregarded.  It  can  be  understood  that  as  the  number  of 
points  at  which  the  boundary  conditions  are  satisfied  approaches 
infinity,  the  actual  boundary  condition  will  be  exactly  satisfied. 
It  is,  however,  possible  to  obtain  a good  approximation  by  choosing 
only  a finite  number  of  points  along  the  edge  x=0  to  impose  the 
constraints[18] [19] [20] . As  will  be  demonstrated  at  the  end 
of  this  section,  five  points  are  usually  good  enough  for  a single 
edge. 

Assume  that,  in  this  case,  I^  points  have  been  chosen  from 
the  edge  x=0  to  impose  the  slope  constraint,  the  following  constraint 
equations  are  obtained: 

Nm  Nn 

a^jm7r/fjsin(n7ry/fy) 

Nm  Nn 

~ ^nin^mn  i ~ 1/2,  ...,Ij  (4.20) 

where  y^  are  the  coordinates  of  the  constraint  points.  The  definition 
Qmni  obtained  directly  from  the  above  equations.  By 
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following  the  same  procedure  demonstrated  in  the  previous  chapter, 
the  Lagrangian  function  is  formed  as  follows: 

„ II 

L = n + S A.G. 

i=l  1 1 


Nm  Nn  o H 

~ ^ ®mn^nin  ” ®mn^mn  ( ^iQmni^mn  ) ] (4.21) 

Taking  partial  derivatives  of  L with  respect  to  each  a^  and  eguating 
the  resulting  eguations  to  zero,  we  obtain  the  following  eguations: 

~ ^mn  “ ^Qmni 

m— 1 , 2 , . . . , Njjj , n=l , 2 , . . . , (4.22) 

and  it  follows  that 


II 


mn 


i?!  ^®inni)  / (^®mn) 


(4.23) 


Substituting  eguation  (4.23)  into  the  constraint  eguations  (4.20)  , 
the  following  set  of  simultaneous  eguations  in  the  Lagrangian 
multipliers  is  obtained: 


II 

2 

i=l 


Nm  Nn 
(22- 
m=l  n=l 


Q Q 

*mn  1 


2E„ 


^mmQm 

m=l  n=l  2E 


i'=l,2, . . . ,Ij  (4.24) 

It  should  be  noted  that  the  number  of  unknowns  are  reduced  from 
Nn,Nn+Ii  to  This  considerably  simplifies  the  solution  of  the 

simultaneous  linear  eguations.  There  are  many  ways  to  solve 
the  eguations,  such  as  the  Gaussian  elimination  method  or  Choleski's 
decomposition  scheme.  The  latter  is  much  more  efficient  than 
the  former,  if  the  system  of  eguations  to  be  solved  is  positive 
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definite[38]  . Besides  the  efficiency,  Choleski  ' s decomposition 
scheme  also  needs  less  computer  memory  than  other  methods  in 
solving  these  kinds  of  problems.  It  can  be  easily  proven  that 
the  system  of  equations  (4.24)  is  positive  definite.  Thus,  Qioleski's 
scheme  is  well  suited  for  their  solution. 

Figure  4.5  shows  a comparison  between  the  theoretical 
deflection  and  three  approximate  deflections  for  a square  plate 
with  a edge  clamped.  The  curves  taken  along  the  section  y=£y/2 
have  been  normalized  by  the  expression  Pq£^Vd.  The  number  of 
terms  for  the  three  approximations  are  11,  15  and  21,  respectively. 
The  number  of  constraint  points  along  the  clamped  edge  is  taken 
to  be  five,  i.e.,  1^=5.  As  can  be  seen  from  figure  4.5,  the 
approximation  curves  approach  the  theoretical  curve  as  the  number 
of  terms  increases.  The  deflection  at  the  centroid  of  the  plate 
is  2 . 65573x10  ^ (Pq£^Vd)  for  the  11-term  approximation,  2. 73920x10'^ (pg^^Vo) 
for  the  15-term  approximation,  and  2 . 79734x10'^  (pg£//D)  for  the 
21— term  ajproximation.  Ihe  difference  between  the  11— term  approximation 
and  the  21-term  approximation  is  approximately  5%  and  that  between 
the  15-term  approximation  and  the  21-term  approximation  is  2%. 
It  can  be  seen  that  the  11-term  approximation  is  good  enough 
in  the  sense  of  convergence.  The  theoretical  deflection  at  the 
centroid  is  0 . 0028  (Pq£^Vd)  [3]  which  is  almost  the  same  as  that 
obtained  from  the  21-term  approximation.  This  shows  that  the 
approximation  is  good  enough  for  practical  application, 
since  the  error  is  only  2%. 


COO-SO'C 
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Figure  4.6  shows  the  deflection  surface  obtained  from  the 
21-term  approximation.  The  deflection  in  this  figure  has  been 
normalized  in  the  same  way  as  above  and  magnified  35  times  relative 
to  the  dimensions  of  the  plate. 


' & qi  - 

Figure  4.6  The  deflection  surface  of  a plate  with  one  edge 
clamped  and  three  edges  simply  supported 

Since  only  discrete  points  along  the  edge  x=0  are  constrained 
to  have  zero  slope,  a variation  of  the  slope  along  this  edge 
between  any  two  successive  constraint  points  can  be  expected. 
This,  however,  is  only  a local  phenomena.  From  figure  4.7,  a 
diagram  plotted  with  the  deflection  at  the  centroid  of  the 
plate  versus  the  number  of  constraint  points,  it  can  seen  that 
as  the  number  of  constraint  points  becomes  greater  than  five, 
the  deflection  of  the  centroid  becomes  nearly  a constant.  As 
was  stated  earlier,  there  will  be  an  exact  clamped  edge  along 
x=0  as  the  number  of  constraint  points  approaches  infinity. 
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Figure  4.7  Deflection  at  the  centroid  versus  number  of 
constraint  points 
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However,  the  deflection  will  be  approximately  the  same  as  the 
case  where  only  five  constraint  points  were  taken.  Obviously, 
the  local  variation  will  not  effect  the  deflection  of  the  plate 
when  a large  enough  number  of  constraint  points  is  taken. 


Plates  with  All  Edges  Clamped 


Using  the  same  strategy,  it  is  also  possible  to  change  the 
previous  case  to  one  with  all  edges  clamped  (figure  4.8).  The 
theoretical  solution  for  this  problem  is  very  tedious  and  it 
is  not  possible  to  have  a general  theoretical  expression  for 
the  deflection[3 ] . However,  the  method  presented  in  this  study 
makes  the  formulation  the  same  as  the  previous  case,  with  the 
exception  of  increasing  the  number  of  constraint  equations. 


Figure  4.8  A clamped  plate 


In  order  to  make  every  edge  clamped,  in  addition  to  equation 
(4.20),  the  following  constraint  equations  must  be  added: 

Nm  Nn 

2^a^n("»V«x)cos(m7r)sin(n7ry/£y)  = 0 


i=Ii+l,i^+2,  . . . ,l2 


(4.25.1) 
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Nm  Nn 

2^a^Jn7r/£y)sin(m7rx/£J  = 0 

i=l2+l,l2+2,  . . . ,13  (4.25.2) 


Nm  Nn 

S^a^Jn7T/£y)  sin(iri7rx/£^)cos(n7r)  = 0 

1=13+1,13+2,  ...  ,14  (4.25.3) 

The  constraint  equation  (4.20)  corresponds  to  the  slope  boundary 
condition  along  edge  x=0  where  Ij  points  are  taken.  Equation 
(4.25.1)  corresponds  to  the  edge  x=£^  where  (I2-I1)  points  are 
taken,  equation  (4.25.2)  corresponds  to  the  edge  y=0  where  (I3- 
I2)  points  are  taken,  and  equation  (4.25.4)  corresponds  to  the 
edge  y=£^  where  (I4-I3)  points  are  taken.  If  is  defined  as 
follows: 


Qmni  = (m7r/£Jsin(n7ry/£y) 

cos  (itiTT)  sin  (nTryy  £y) 
(nTT/fy)  sin(m7rx/£^) 

(n7T/£y)  sin  (m7rx/€^)  cos  (nvr) 


if  i < Ij 
if  li  < i < I2 
if  I2  < i < I3 
if  I3  < i < I4 
(4.26) 


then  all  formulations  will  be  the  same  as  the  previous  case, 
except  that  the  number  of  constraint  equations  is  I4. 

Since  it  is  almost  impossible  to  obtain  a general  theoretical 
deflection  formulation  for  this  problem,  the  result  obtained 
from  the  finite  element  method  will  be  taken  as  a reference  for 
comparison.  The  finite  element  method  is  commonly  considered 
an  accurate  numerical  scheme  and  is  often  used  as  a comparison 
reference [17]  [39]  . For  a clamped  square  plate,  the  theoretical 
maximum  deflection,  which  occurs  at  the  centroid,  can  be  computed 
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to  be  0. 00126  (Pq£^Vd)  and  that  for  the  finite  element  method  is 
0. 00127  (Pq£^Vd)  . The  difference  is  less  than  1%,  and  shows  the 
accuracy  of  the  finite  element  method  in  this  type  of  problem. 
It  should  be  noted  that  both  the  finite  element  method  and  the 
modified  Rayleigh-Ritz  method  provide  only  approximate  solutions. 
In  general,  these  approximations  are  not  the  same  but  will  be 
close  to  one  another. 

Taking  a square  plate  as  an  example,  figure  4.9  shows  a 
comparison  between  the  deflections  obtained  from  the  finite  element 
method  with  400  square  elements  and  three  approximations  using 
and  21  terms  respectively.  The  deflection  curves 
shown  in  the  figure  are  taken  along  the  section  y=£y/2  and  have 
been  normalized  by  the  expression  Pq£^Vd.  As  can  be  seen  from 
the  plot,  convergence  occurs  when  the  number  of  terms  reaches 
21.  The  difference  between  the  15—  and  21— term  approximations 
is  about  3.6%,  and  the  difference  between  the  21-term  approximation 
and  the  result  obtained  from  the  finite  element  method  is  about 
0.5%.  The  21-term  approximation  can  be  regarded  as  a very  good 
approximation  for  engineering  applications. 

Figure  4.10  shows  the  deflection  surface  obtained  from  the 
21-term  approximation.  The  deflection  in  the  figure  has  also 
been  normalized  in  the  same  way  as  above  and  magnified  80  times 
with  respect  to  the  dimensions  of  the  pate. 
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Figure  4.9  Deflection  curves  of  a clamped  plate  taken 
along  the  section  y=£  /2 
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Figure  4.10  The  deflection  surface  of  a clamped  plate 

Nonrectanqular  Plates 

One  of  the  advantages  of  the  proposed  method  is  that  it 
is  possible  to  obtain  the  deflection  surface  of  a nonrectangular 
plate  from  the  formulation  of  a rectangular  plate.  The  concept 
of  approximating  a nonrectangular  plate  from  a rectangular 
plate  requires  imposing  suitable  constraints  to  the  interior 
of  the  rectangular  plate  so  as  to  make  the  geometrical  situation 
at  the  constrained  points  the  same  as  what  is  to  be  expected 
along  the  boundary  of  the  nonrectangular  plate. 

Taking  a circular  plate  as  an  example  (shown  in  figure  4.11) , 
a circle  is  "drawn"  inside  the  rectangular  plate  (a  square  in 
this  case)  and  constraint  points  are  taken  along  the  circle. 
The  two  constraints  which  need  to  be  imposed  on  the  chosen  points 
are  that  the  deflection  at  these  points  should  be  zero,  and  that 
the  slope  in  the  normal  direction  to  the  circle  should  be  zero. 
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Figure  4.11  A clamped  circular  plate 


If  the  same  series  as  equation  (4.12)  is  used,  and  the  applied 
load  is  uniform,  the  above-mentioned  constraints,  in  mathematical 
form,  are 


G,  = 


Nm 

v(x^,y.)  = S 
' m=l 


Nn 

^mnSinCmTTx/fJ  sin(n7ry/fy) 


and 


i=l,2. 


Ii.  j=i 


(4.27) 


Gj  - 3v/5n  = cos(^)3v/3x  + sin(^)3v/3y 


Nm  Nn  cos(0)m7T  mTTX.  nny 

= „?,  COS(— )sin(— -^) 

^ Y V ^ 


m=l  n=l 


sin(5)n7T  mTTX.  nTry. 

+ sin( ^)cos( ^)] 

i f P 


1 l,2,...,Ij^,  j — 1+Ij^  (4.28) 

where  6 is  the  angle  between  the  outward  normal  to  the  circle 
and  the  x axis.  In  equations  (4.27)  and  (4.28),  x^  and  y^  are 
the  coordinates  of  the  constraint  points  and  I^  is  the  number 
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of  constraint  points.  Since  two  constraints  are  imposed  on  each 
constraint  point,  the  total  number  of  constraint  equations  is 
21^. 

If  Q^rii  defined  as  follows: 

Qmni  = sin(m7Tx/£Jsin(n7ry/£y)  if  i < Ij 

cos(0)m7r  mTTx.  mry 

= cos( ^)sin( 

i Z P 

sin(0)n7r  myrx-  n7ry 

+ sin(— - — )cos( -)  if  i > I, 

^ ^ ' 

(4.29) 

then  the  formulation  will  take  the  same  form  as  in  the  previous 
cases  and  the  deflection  surface  can  be  obtained  in  the  same 
manner. 

Taking  a square  plate  as  the  original  plate  and  imposing 
constraints  on  a circle  whose  center  is  the  centroid  of  the  square 
and  whose  radius  is  equal  to  half  of  the  length  of  the  plate 
(figure  4.11) , the  result  shown  in  figure  4.12  can  be  obtained. 
The  deflection  curves  shown  in  figure  4.12  have  been  normalized 
by  the  expression  These  curves,  along  with  the  theoretical 

curve[3] , are  taken  along  a diameter  of  the  circle.  The  numbers 
of  terms  for  the  three  approximations  are  15,  17,  and  21,  respectively. 
The  difference  of  the  maximum  deflection  between  the  15  and  21 
term  approximations  is  7.3%  but  that  between  the  17  and  21  term 
approximations  is  only  1.5%.  The  convergence  can  be  considered 
to  occur  when  the  number  of  terms  is  21.  As  compared  to  the 
theoretical  value,  the  error  of  the  maximum  deflection  of  the 
21-term  approximation  is  less  than  3%. 
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Figure  4.12  Deflection  curves  of  a clamped  circular  plate 
taken  along  a diameter 
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In  figure  4.13,  the  deflection  surface  obtained  with  the 
21-term  approximation  is  shown.  The  deflection  in  the  figure 
has  been  normalized  in  the  same  way  as  above  and  magnified  100 
times  with  respect  to  the  dimensions  of  the  plate. 


* CX)  * 

Figure  4.13  The  deflection  surface  of  a circular  plate 
Plates  with  Mixed  Boundary  Conditions 

The  problem  of  plates  with  mixed  boundary  conditions  is 
a very  difficult  one  in  the  classical  theory  of  plates.  Generally, 
the  finite  element  method  and  the  finite  difference  method  are 
applied  to  solve  this  type  of  problem.  Since  the  proposed  method 
allows  the  boundary  conditions  to  be  satisfied  through  the  use 
of  Lagrangian  multipliers  on  a finite  number  of  discrete  points 
along  the  boundary,  it  is  possible  to  handle  a mixed  boundary 
condition  problem  without  any  major  modification  of  the  established 
formulations . 
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The  strategy  of  applying  the  proposed  method  to  a mixed 
boundary  condition  problem  is  to  impose  suitable  constraints 
according  to  the  desired  boundary  conditions  at  every  constraint 
point  where  the  boundary  conditions  are  not  completely  satisfied 
by  the  approximation  series,  and  to  disregard  those  points  where 
boundary  conditions  are  satisfied.  The  situation  shown  in  figure 
4.14  is  taken  as  an  example  for  illustration. 
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Figure  4.14  A plate  with  mixed  boundary  conditions 


The  constraint  points  are  chosen  inside  the  following 
regions  along  the  edges  of  the  plate: 
region  1:  x=0  and  0 < y < 
region  2:  x=f^  and  fy/2  < y < 
region  3:  y=0  and  0 < x < Z^/4 
region  4:  y=fy  and  f^2  < x < 

With  the  same  approximation  series  as  equation  (4.12) , the  constraint 
imposed  on  regions  1 and  2 is  dv/dx=0,  and  that  on  regions  3 
and  4 is  5v/3y=0.  With  a similar  definition  of  as  equation 
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(4.26),  the  formulations  will  be  the  same  as  in  the  previous 
cases.  The  general  procedure  is  also  the  same. 

Figure  4.15  shows  a comparison  between  the  result  obtained 
from  the  finite  element  method  with  400  elements  and  two  approximations 
with  the  number  of  terms  egual  to  21  and  23.  Since  it  is  impossible 
to  obtain  a theoretical  solution  for  this  case,  the  result  of 
the  finite  element  method  is  taken  as  a reference  with  the  same 
reason  as  that  stated  previously.  In  this  example,  the  length 
of  the  plate  is  taken  to  be  twice  the  width  The  deflection 
curves,  which  are  taken  along  the  section  x=£^/2,  have  been  normalized 
by  the  expression  Pq^^^Vd.  The  difference  between  the  21-  and 
23-term  approximations  is  about  2.2%  and  that  between  the  23- 
term  approximation  and  the  finite  element  method  is  4.7%.  It 
can  be  seen  from  the  figure  that  the  left  side  is  simply  supported 
and  the  right  side  is  clamped,  as  reguired  by  figure  4.14.  This 
can  also  be  seen  from  figure  4.16,  which  shows  the  entire  deflection 
surface  obtained  with  the  23-term  approximation.  Figure  4.16 
has  also  been  normalized  in  the  same  way  as  figure  4.15  has. 

In  figure  4.16,  the  trace  of  the  simple  support  along  the 
edges  parallel  to  the  x axis  can  be  easily  identified,  but  is 
not  readily  seen  on  the  edge  y=£y.  Since  £y  is  only  half  of  £^, 
the  edges  parallel  to  the  y axis  resist  deflection  more  than 
the  edges  parallel  to  the  x axis. 
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Figure  4.16  The  deflection  surface  of 
a mixed  boundary  plate 

Stiffened  Plates 

The  problem  of  a stiffened  plate  can  also  be  conveniently 
handled  by  the  modified  Rayleigh-Ritz  method.  The  stiffeners, 
which  are  in  the  form  of  beams  in  this  investigation,  store  energy 
when  they  are  bent.  The  bending  of  the  stiffeners  is  the  same 
as  the  beam  bending  which  has  been  discussed  in  the  previous 
chapter.  The  deflections  of  the  stiffeners  should  be  compatible 
with  the  deflection  of  the  plate.  This  constitutes  the  constraints 
between  the  plate  and  the  stiffeners. 

As  a demonstration  of  how  to  apply  the  proposed  method  to 
a stiffened  plate,  a simply-supported  square  plate  reinforced 
by  only  one  stiffener  is  considered,  as  shown  in  figure  4.17. 

If  the  deflection  function  of  the  stiffener  is  taken  as 
the  equation 
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^ Nk  k7T^ 

= S^dkSin(— 


(4.30) 


where  is  the  length  of  the  stiffener  and  ^ is  the  local  coordinate 
describing  the  stiffener,  then  the  strain  energy  stored  in  the 
stiffener  is 


U 


S 


EItt"* 


Nk 


2i.4 


■S  d/k 
k=l 


(4.31) 


which  is  similar  to  equation  (3.6).  The  strain  energy  of  the 
stiffener  should  be  added  to  the  Lagrangian. 


Figure  4.17  A stiffened  plate 

The  relationships  between  the  coordinate  ^ and  x and  y are 
^cos(e)  = X - Xq  (4.32.1) 

^sin{9)  = y - yg  (4.32.2) 

where  Xq  and  yg  are  the  coordinates  of  the  origin  of  the  local 
coordinate  £ in  terms  of  the  coordinate  system  x and  y.  The 
deflection  compatibility  requires  that  v(x,y)  = u(0  along  the 
stiffener.  If  the  deflection  surface  of  the  plate  is  approximated 
by  the  same  sine  series  as  equation  (4.12)  , then  choosing  a finite 
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number  (I^)  of  discrete  points  along  the  stiffener  will  result 
in  the  following  constraint  equations: 

Gi(a„n'd^)  = v(x^,y^)  - u(e^)  = 0 

i=l,2,..,Ii  (4.33) 

where  x^,  and  follow  the  relationships  of  (4.32).  These 
equations  can  be  rewritten  as 


Nm  Nn 

S 2 a ( 

m=l  n=l 


Nk 

- = 


i=l,2,...,Ii  (4.34) 


where 


and 


= sin(k7rCi/«^) 


Qmni  = sin(m7rx/£^)sin(n7ry/fy) 
The  Lagrangian  L is 


(4.35.1) 


(4.35.2) 


L = n + S 

i=l  ’ ’ 


Nm  Nn 

= 22  [E, 

m=l  n=l  ' 


mn^mn 


- S a 

mn  mn 


II 


Nk 


II 


+ 2(A  Q a ) 1 - 2 du  ( 2 X T^  ) 

i^mm  mn'  J k=l  ^ 'iii  i ki  ' 


(4.36) 

Following  the  same  procedure  as  the  previous  cases,  and  taking 
partial  derivatives  of  the  Lagrangian  function  L with  respect 
to  all  parameters  a^^  and  d|^,  the  following  relationships  between 
the  parameters  and  the  Lagrangian  multipliers  are  obtained: 


^mn 


2E„ 


- 2 

i=l  ^ 


2E„ 


(4.37) 


and 
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d 


k 


II 

2 


k=l 


[ Th  / (2F,)  ] 


(4.38) 


where 


4 


(4.39) 


Substituting  equations  (4.37)  and  (4.38)  into  equation  (4.34), 
the  following  equations  in  are  obtained: 


^ ^ gmA  l- 

i=l  ' ni=l  n=l  2E 


+ E 
k=l 


Nk  T^jTkj. 


2F, 


Nm  Nn  S„Q_ 


2 2 
in=l  n=l  2E 


i'  l,2,...,Ji 


(4.40) 


Solving  for  the  undetermined  multipliers  in  the  above  equations 
and  substituting  them  into  equations  (4.37)  and  (4.38),  the  parameters 
a^^  and  d^  are  found.  With  these  parameters,  the  deflection  of 
the  plate  can  be  readily  calculated. 

Since  it  takes  more  terms  for  the  proposed  method  to  meet 
the  requirement  of  unsatisfied  boundary  conditions  (by  the  assumed 
deflection  function)  than  to  satisfy  the  compatibility  requirement 
between  the  stiffener  and  the  plate,  less  terms  are  necessary 
to  give  a convergent  solution  for  the  present  case  (where  the 
boundary  conditions  are  all  simply  supported  and  are  automatically 
satisfied  by  the  assumed  deflection  function)  than  for  some  of 
the  previous  cases  with  boundary  conditions  not  all  satisfied 
by  the  assumed  deflection  function,  for  example,  the  plate  with 
mixed  boundary  conditions.  Therefore,  it  takes  only  nine  terms 
to  obtain  a convergent  result  for  the  case  of  a square  plate 
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reinforced  by  a stiffener  with  Xq=(.^/2  and  yo=0,  and  6 equal 
to  90  degrees.  Figure  4.18  shows  the  deflection  surface  of  the 
stiffened  square  plate.  The  relative  rigidity,  which  is  defined 
as  the  product  of  the  flexural  rigidity  of  the  plate  times  the 
thickness  of  the  plate  divided  by  the  bending  rigidity  of  the 
stiffener  (Dh/EIg)  is  computed  to  be  e.lsexlO'"^.  In  figure  4.18, 
the  trace  of  the  stiffener  can  be  easily  identified.  The  deflection 
surface  has  been  normalized  in  the  same  manner  as  the  previous 
cases. 


Figure  4.18  The  deflection  surface  of  a plate  reinforced 
by  a stiffener  parallel  to  the  y axis 

If  there  is  more  than  one  stiffener,  the  strain  energy  of 
each  stiffener  should  be  added  to  the  total  potential  energy. 
Each  stiffener  should  be  assumed  to  have  a different  deflection 
function  with  a different  set  of  unknown  displacement  parameters, 
and  the  compatibility  between  the  plate  and  each  stiffener 
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should  be  regarded  as  constraints.  The  number  of  sets  of  parameters 
will  be  equal  to  the  number  of  stiffeners,  plus  one  set  of  parcuneters 
for  the  plate  itself.  Relationships  between  each  set  of  parameters 
and  the  undetermined  multipliers  are  obtained  by  differentiating 
the  Lagrangian  with  respect  to  each  parameter  of  each  set.  Combining 
these  relationships  with  the  constraints  of  deflection  compatibility, 
a set  of  simultaneous  equations  similar  to  equation  (4.40)  is 
obtained.  Thus  the  deflection  of  the  plate  can  be  obtained. 

Figure  4.19  shows  the  deflection  surface  of  a square 
plate  diagonally  reinforced  by  two  stiffeners.  For  this  case, 
a convergent  result  can  be  obtained  with  a nine-term  approximation. 
The  relative  rigidity  is  l.OxlO"'*  for  both  stiffeners.  The  traces 
of  the  stiffeners  are  easy  to  identify  from  the  figure.  The 
deflection  surface  has  been  normalized  as  before. 


Figure  4.19  The  deflection  surface  of  a 
diagonally  stiffened  plate 
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Although  the  stiffened  plate  exctmples  shown  above  are  simple, 
the  application  of  the  proposed  method  to  this  type  of  problem 
can  clearly  be  extended  to  more  complex  geometries.  In  the  following 
section,  a demonstration  of  how  to  handle  a stiffened  double 
plate,  called  a honeycomb  sandwich  plate,  will  be  considered 
to  illustrate  the  versatility  of  the  proposed  method. 

Honeycomb  Sandwich  Plates 

It  is  extremely  difficult  for  conventional  methods,  even 
for  the  finite  element  method,  to  handle  the  problem  of  honeycomb 
sandwich  plates.  Neither  the  theoretical  solution  nor  a finite 
element  method  solution  is  therefore  available  for  the  case  shown 
in  figure  4.20.  However,  by  extending  the  procedure  used  to 
analyze  the  stiffened  plates  discussed  in  the  previous  section, 
it  is  possible  to  conveniently  deal  with  this  type  of  problem 
with  the  proposed  method.  This  is  one  of  the  advantages  of 
the  proposed  method  over  conventional  methods. 


Figure  4.20  A honeycomb  sandwich  plate 
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A sandwich  plate  can  be  regarded  as  a combination  of  three 
parts,  two  of  them  are  the  top  and  bottom  plates  and  the  third 
represents  the  stiffeners.  Since  it  was  pointed  out  at  the  outset 
of  this  chapter  that  the  plate  problems  dealt  with  in  this  study 
are  considered  linearly  elastic,  the  principle  of  superposition 
is  applicable.  The  deflection  of  the  top  or  bottom  plate  of 
the  sandwich  structure  can  then  be  viewed  as  a superposition 
of  two  parts.  One  is  the  overall  deflection  of  the  whole  structure 
and  the  other  is  the  deflection  of  the  top  or  bottom  plate  when 
the  stiffeners  are  assumed  rigid. 

The  deflecticn  of  either  plate  can  then  be  obtained  by  stperiirposing 
their  "local"  deflection  on  the  "global"  deflection.  The  terminologies 
"local"  and  "global"  are  used  as  an  indication  of  the  deflection 
of  the  plates  when  the  stiffeners  are  assumed  rigid,  and  the 
overall  deflection  of  the  whole  structure.  Kirchhoff's 
assumptions  for  the  analysis  of  plates  will  be  taken  to  be  valid. 
The  strain  energy  stored  in  both  the  top  plate  and  the  bottom 
plate  thus  consists  of  three  parts.  The  first  is  due  to  the 
local  deflection  when  the  global  deflection  is  zero,  the  second 
is  due  to  the  global  deflection  when  the  local  deflection  does 
not  exist,  and  the  third  is  the  interaction  between  the  local 
and  global  deflections.  Although  we  are  dealing  with  a linear 
system,  it  is  not  possible  to  apply  the  principle  of  superposition 
to  energy.  It  is  for  this  reason  that  the  third  part  of  the 
strain  energy  is  necessary.  The  third  part  can  be  regarded  as 
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a gain  in  strain  energy  when  the  stress  produced  by  the  local 
deflection  undergoes  the  global  deflection. 

To  treat  the  honeycomb  structure  as  three  parts,  the  deflection 
functions  for  the  top  plate,  the  bottom  plate,  and  the  whole 
structure  are  independently  assumed.  The  overall  deflection 
function  of  the  whole  structure  can  be  regarded  as  the  deflection 
of  a pseudo  plate  which  is  located  on  the  middle  plane  of  the 
whole  structure,  and  which  is  parallel  to  both  the  top  and  bottom 
plates.  The  deflection  functions  for  the  top  and  bottom  plates 
cover  only  their  local  deflections.  Therefore  the  deflection 
function  of  the  top  plate  on  its  intersection  with  the  stiffeners 
has  a zero  value,  and  the  same  is  true  for  the  bottom  plate. 
This  is  the  constraint  for  the  local  deflection  functions  of 
both  the  top  and  the  bottom  plate.  The  real  deflection  of  the 
top  or  bottom  plate  is  the  sum  of  its  local  deflection  and  the 
deflection  of  the  pseudo  plate. 

To  follow  the  same  rule  of  calculating  the  potential  energy 
and  forming  the  constraint  eguations  as  was  used  in  the  previous 
cases,  the  local  deflection  functions  for  the  top  and  bottom 
plates  are  taken  respectively  to  be 

Nm  Nn 

Vt(x,y)  b^„sin(m7rx/€Jsin(n7ry/£y)  (4.41) 

and 

Nm  Nn 

Vb(x,y)  =^S  c^^sin(m7Tx/£Jsin(n7ry/£y) 


(4.42) 
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The  subscripts  t and  b denote,  respectively,  the  top  plate  and 
the  bottom  plate.  The  global  deflection  function  is  taken  to 
be  the  same  as  equation  (4.12).  The  deflection  functions  of 
the  stiffeners  are  taken  as 

> N n7r£. 

^ i=l,2,...,I^  (4.43) 

^ i 

where  1^.  is  the  number  of  stiffeners,  is  the  local  coordinate 
describing  the  ith  stiffener,  and  is  the  length  of  the  ith 
stiffener. 

The  strain  energy  due  to  the  local  deflection  of  the  top 
plate  is 


Nm  Nn  , mTT  , n7T 

Um  = — s b ^ r( )^  + ( 

8 m=l  n=l  ^ ^ ' ' £ 

X y 


and  that  of  the  bottom  plate  is 

Nm  Nn  , mTT 


(4.44) 


Ubl  = 


mil  mi  o lU/l  o 

s c ^ r( + ( 

8 m=l  n=l  ^ '•  ' £ ^ ^ 


nTT 


(4.45) 


The  strain  energy  due  to  the  global  deflection  of  the  top  plate 
is 


U,  = 


D £ £ Nm  Nn  p 
^ ^ ^ i:  S a ^ 

8 m=l  n=l  ™ 


[( 


mTT 


)'  + 


mr 


(4.46) 


where  is  defined  as 

Eh  (H-h)^ 

Dc  = D + 2 (4.47) 

(1-1^  ) 4 

The  symbol  H denotes  the  overall  height  of  the  honeycomb  structure. 
Because  of  the  symmetry  about  the  middle  plane  of  the  v^ole  structure, 
the  strain  energy  due  to  the  global  deflection  of  the  bottom 
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plate  is  the  same  as  equation  (4.46).  The  third  part  of  the 
strain  energy  for  the  top  plate  is 


Ut3  = 


Nm  Nn 


mTT 


2 2 a„„b„„  [ ( ) ^ + ( 

4 m=l  n=l  ^ £ ' ' 


n7T 


(4.48) 


and  that  for  the  bottom  plate  is 


Df  f Nm  Nn  iu«  5 

Ub3  = - ' ' 


2 2 a_-C_„  [ ( ) + ( 

4 m=l  n=i  mn  mn  L k / V 


mTT 


n7T 


(4.49) 


The  strain  energy  of  each  stiffener  is 


Usi  = 


EItt^  Nk 


4f, 


^2  d 

3 '^ik 


i 1,2, ...,Xp 


(4.50) 


i k=l 


If  the  external  load  applied  on  the  top  surface  of  the  structure 
is  uniformly  distributed  with  intensity  p^.,  then  the  work  done 
by  this  load  is 

Pf^x'^y  ^ 


Wt  = 


1-cos  (mTT)  1-cos  (nTT) 
■^,11  (amn+b.J  ( ) ( ^ ) 


m 


n 


(4.51) 

Similcirly,  the  work  done  by  a uniform  load  with  intensity  P|^  distributed 
over  the  bottom  surface  of  the  structure  is 


Wb  = 


Pb'^x'^v 


m=l  n=l 


2 2 (a„„+c„„)  ( 

,_i  -._i  ' mn  mn'  ' 


l-cos(mTr)  1-cos  (nTT) 


m 


) ( 


n 


(4.52) 


Therefore,  the  total  potential  energy  is 


Ir 


n = + 2u,  + Ut,  + + 2 u.,  - w,  - w, 


t3  ^ '-’b3 


i=l 


(4.53) 
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The  constraints  between  the  global  deflection  and  the  deflections 
of  the  stiffeners  are 


Giij  = 


Nm  Nn 
S 2 
m=l  n=l 


Nk 

“ k?i  '^ik'^kij  = 


i 1,2,...,!^,  j— 1,2,...,J^  (4.54.1) 


The  constraints  for  the  deflection  function  of  the  top  plate 
are 


Nm  Nn 

G,..  = 2 2 b Q =0 

n]=l  n=2  mn i j 

i=l,  2 

and  for  the  bottom  plate  are 


1^,  j=l,2,...,J.  (4.54.2) 


Nm  Nn 

G,..  = 22cQ  . =0 

m=l  n=l 

l~l/2,...,I^,  j~lf2,...,J^  (4.54.3) 

where 

Qmnij  = sin(m7rx,j/£^)sin(n7ry,j/£y)  (4.55) 

and 

Tkij  = sin(k7re/£.)  (4.56) 

The  symbol  J.  denotes  the  number  of  constraint  points  taken  on 
the  ith  stiffener. 

The  Lagrangian  function  takes  the  form: 


L = n 


Ir  Ji 

+ 2 2 ( 

i=l  j=l  ' 


+ + f..G3,.  ) 


(4.57) 


By  taking  partial  derivatives  of  L with  respect  to  all  displacement 
parameters  a^^,  b^^,  c^^,  and  d^,^  (i=l,  2 , . . . , I^)  , and  substituting 
the  results  into  the  constraint  equations  (4.54)  , the  following 
equations  are  obtained: 
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Ir  Ji  Nm  Nn  QmninQmni'i’ 

{ S S_,[  ( .5  g_7  ^ ) ] } 


(D^-D)£^£y  i=l  j=l  ’J  m=l 

2^1. ^ Ji-  ^ Nk 

V«  ' )i  = ° 


EItt 


-2 


i'  1/2,...,!^,  3'~1/2,...,J^ 


Ir  Ji  ^ Nm  Nn  QmniiQmni'i’ 

S E,  _s  _s — 


(4.58.1) 


(D^,-D)£^£^  i=l  j=l  n=l  m=l  £, 

+ [ 


-x'-y 

2 


4 Ir  Ji  Nm  Nn  QmniiQmni'i' 

+ ] S,  S,  (Mii  _S  2; nmi_JD!lLJ_  ) 


(Dj,-D)«^£y  i=l  j=l  n=l  m=l 


'x'y 


S 2,  (f  ii  iT  _s-  ) 


mn 


(D^-D)£^£y  i=l  j=l  n=l  m=l 


4Pt  ( 1-COS  (mTT)  ) ( 1-COS  (n7T)  ) 

Dtt^  m n E 


i'  1/2,...,!^,  j'— 1,2,...,J^. 


(4.58.2) 


-2 


Ir  Ji  . Nm  Nn  Qmni  -iQmni  ' i ' 

2,  S,  (A..  _S  S; Eml^mn^^j 


(D^-D)  £^€y  i=l  j=l  n=l  m=l 


+ 


Ir  Ji  Nm  Nn  QmniiQmni'i' 

2^  (/i.  . _S  S; nmLLJlILJ— ) 


+ [ 


(D^-D)£^£y  i=l  j=l  n=l  m=l 

2 


mn 

4 Ir  Ji  Nm  Nn  QmniiQi 


(D^-D)£^£y  ^ i?if=i^^ij  n^im=l 


mn  1 .1  ^mn  i .r 


'x^y 


4Pb  ( 1-cos  (lUTT)  ) (1-cos  (n7T)  ) Qn,ni-i- 

Dtt^  III  n E„„ 

mn 

i'=l,2,  . . . ,1^,  j '=1,2,  . . . ,J^.  (4.58.3) 

Solving  equation  (4.58)  gives  the  undetermined  multipliers 
A^j,  ju.j,  and  With  these  multipliers,  the  displacement  parameters 
can  be  calculated  by  taking  partial  derivatives  of  equation  (4.57) 
with  respect  to  each  parameter.  The  computer  program  used  to 
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solve  the  system  of  equations  in  A^j,  and  and  to  generate 
the  deflection  surface  is  listed  in  the  appendix  for  reference. 

Figures  4.21  and  4.22  show  the  deflection  surfaces  of  a 
sandwich  plate  with  a length  of  3 units,  a width  of  1 unit,  and 
an  overall  height  of  0.1  units.  The  thickness  of  the  plates 
is  0.01  units  and  that  of  the  stiffeners  is  0.01  units.  Figure 
4.21  shows  the  deflection  surface  of  the  top  plate.  The  deflection 
in  this  figure  is  magnified  3,000  times  with  respect  to  the  dimensions 
of  the  structure.  Figure  4.22  shows  the  deflection  surface  of 
the  bottom  plate.  The  deflection  of  this  figure  is  magnified 
30,000  times  with  respect  to  the  dimensions  of  the  structure. 
It  can  be  seen  that  the  overall  deflection  (figure  4.22)  is  much 
smaller  than  the  local  deflection  of  the  top  plate  (figure  4.21)  . 

In  figure  4.23,  the  maximum  deflections  of  the  top  plate 
the  pseudo  plate  are  shown.  As  expected,  the  difference 
between  these  two  curves  decreases  as  the  thickness  of  the  plate 
increases  and  approaches  half  of  the  overall  height  of  the  honeycomb. 
The  load  in  this  case  is  1 force  unit.  Poisson's  ratio  is 
taken  to  be  0.3,  and  Young's  modulus  is  30x10®  force  unit  over 
length  unit  square. 
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Figure  4.21  The  deflection  surface  of  the  top  plate 


Figure  4.22  The  deflection  surface  of  the  bottom  plate 
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CHAPTER  5 

APPLICATIONS  TO  PLATES  WITH  HOLES 
Introduction 

A plate  with  arbitrarily  shaped  holes  is  a coitiition 
structural  element,  but  is  extremely  difficult  to  handle  by 
the  classical  theory  of  plates.  The  finite  element  method, 
so  far,  is  the  approach  that  analysts  prefer  to  use.  However, 
as  pointed  out  previously,  the  disadvantages  of  the 
preprocessing  and  postprocessing  effort  makes  the  finite 
element  method  an  expensive  way  to  deal  with  simple  structural 
elements. 

The  Rayleigh-Ritz  method  using  undetermined  multipliers 
can  be  effectively  applied  to  this  type  of  plate  problem.  In 
the  example  of  nonrectangular  plates  discussed  in  Chapter  4 , 
the  clamped  boundary  is  moved  inward  from  a rectangular  plate 
to  match  the  boundary  curve  of  the  nonrectangular  plate,  which 
is  artificially  drawn  inside  the  rectangular  plate.  This 
concept  can  be  further  extended  to  the  problem  of  plates  with 
holes. 

Taking  a hole  with  a clamped  boundary  as  an  example,  a 
pseudo  plate  of  the  same  size  as  the  hole  can  be  fitted  into 
the  hole,  with  its  boundary  clamped.  This  pseudo  plate  can 
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then  be  integrated  with  the  original  plate,  inside  which  the 
hole  is  located,  and  analyzed  as  a whole  with  the  same 
procedure  shown  in  the  previous  chapter.  The  pseudo  plate 
inside  the  hole  can  then  be  disregarded.  Under  certain 
situations,  such  as  a hole  with  free  edges,  it  is  necessary 
to  superimpose  additional  forces  (generalized  forces)  along 
the  boundary  so  as  to  comply  with  the  force  boundary 
conditions.  It  has  been  shown [40]  that  load  sources  can  be 
applied  externally  to  the  plate  boundary  to  satisfy  the 
boundary  conditions  at  a finite  number  of  points  along  the 
boundary.  This  concept  is  incorporated  to  solve  the  type  of 
problem  to  be  discussed  in  the  following  sections.  With  this 
concept,  forces  are  artificially  superimposed  inside  the  hole 
but  very  close  to  its  boundary,  so  as  to  make  the  boundary  act 
as  desired.  As  stated  in  Chapter  4,  a variation  between  any 
two  successive  constraint  points  is  inevitable  but  this 
variation  is  only  a local  phenomena.  The  following 
investigation  is  also  based  on  Kirchhoff's  three  assumptions 
as  outlined  in  Chapter  4. 

An  illustration  of  how  to  handle  the  problem  of  a square 
plate  with  a square  clamped  hole  will  be  given  first.  This 
is  a simple  problem  because  the  clamped  edges  have  already 
provided  the  necessary  forces  to  satisfy  the  boundary 
conditions  and  no  additional  forces  are  necessary. 
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Scmare  Plates  with  Clamped  Square  Holes 

The  problem  is  described  in  figure  5.1.  The  edges  of 
the  hole  are  clamped  and  the  boundary  of  the  plate  is  simply 
supported . 


y 


s.s, 


Figure  5.1  A plate  with  a clamped  hole 


As  stated  in  the  previous  section,  a pseudo  plate  is 
fitted  into  the  hole  and  integrated  with  the  original  plate. 
The  plates  are  then  analyzed  as  a whole  with  the  same 
procedure  as  that  discussed  in  Chapter  4 . By  taking  the  same 
approximation  series  as  equation  (4.12),  the  formulation  of 
the  strain  energy  of  the  plate,  including  the  pseudo  plate 
inside  the  hole,  is  the  same  as  equation  (4.13) . If  the  load 
acting  upon  the  plate  is  uniform  over  the  whole  surface, 
including  the  pseudo  plate,  with  the  intensity  Pq,  then  the 
formulation  of  the  work  done  by  the  external  load  is  the  same 
as  equation  (4.14).  The  constraints  are  zero  deflection  and 
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zero  normal  slope,  i.e.,  v=0  and  3v/3n=0,  along  the  edges  of 
the  hole.  The  formulation  will  be  the  same  as  that  used  in 
the  case  of  nonrectangular  plates  provided  that  is  defined 
as  eguation  (4.29). 

For  the  proposed  method,  a 20-term  approximation  gives 
a convergent  result.  Figures  5.2  and  5.3  show  comparisons  of 
the  result  obtained  from  the  20-term  approximation  with  that 
provided  by  the  finite  element  method  with  300  elements.  Both 
figures  have  been  normalized  by  the  expression  Figure 
5.2  is  taken  along  the  section  y=£y/2.  Figure  5.3  is  taken 
along  the  section  y=0.15£y.  The  maximum  deflection  occurs 
near  x=0.15£jj  and  y=0.15£y.  The  maximum  deflection  obtained 
from  the  finite  element  method  is  3 . 10725x10'^  (Pq£^Vd)  , and 
that  obtained  from  the  proposed  method  is  2 . 96271x10'®  (P(j£//D)  . 
The  difference  is  less  than  5%. 

Figure  5.4  shows  the  deflection  surface.  The  clamped 
edges  of  the  hole  can  be  identified  from  the  plot.  The 
deflection  in  figure  5.4  has  been  normalized  in  the  same 
manner  as  figures  5.2  and  5.3  and  magnified  1200  times  with 
respect  to  the  dimensions  of  the  plate. 

Simulation  of  a Cantilever  Beam  From  a Clamped  Beam 

Although  it  takes  no  extra  effort  to  handle  a clamped 
hole,  a more  sophisticated  technique  is  necessary  to  handle 
a hole  with  free  edges.  Before  solving  a plate  with  a free 
edge  hole,  an  analogous  beam  problem  will  be  considered. 
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Figure  5.4  The  deflection  surface  of  a plate 
with  a clamped  hole 

consisting  of  the  simulation  of  a cantilever  beam  from  a 
clamped-clamped  beam.  Consider  the  following  condition 
(figure  5.5)  where  the  middle  portion  of  the  beam  is  not 
loaded,  and  the  other  portions  of  the  beam  are  uniformly 
loaded  with  load  intensity  Pg. 


V 


Po 

Po 

% 

f 

1 

X 


Figure  5.5  A clamped-clamped  beam  partially  loaded 


By  taking  the  same  series  as  eguation  (3.13)  to  be  the 
approximate  deflection  function,  the  strain  energy  will  take 
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the  same  form  as  equation  (3.10)  with  the  upper  limit  of  the 
summation  replaced  by  N.  The  work  done  by  the  applied  load 
is 


N Pof 

W = S a ( 1-cos  (nTT/4 ) +COS  (3n7T/4 ) -cos  (nTT/4 ) ) 

n=i  nir 

(5.1) 


Following  the  same  procedure  as  illustrated  in  Chapter  3,  the 
deflection  curve,  which  is  labeled  iteration  1 in  figure  5.6, 
can  be  obtained.  The  bending  moment  of  the  beam  is  shown  in 
figure  5.7,  also  labeled  iteration  1.  In  the  same  plots, 
another  curve  labeled  "double  cantilever"  is  drawn 
corresponding  to  the  beam  shown  in  figure  5.5,  except  that  it 
is  cut  off  at  the  center  to  simulate  two  cantilever  beams. 
This  configuration  will  be  referred  to  as  the  "double 
cantilever  beam."  Note  that  when  the  clamped-clamped  beam  is 
cut  off  at  the  center,  both  halves  can  be  treated  as  a 
uniformly  loaded  cantilever  beam  of  length  £/4  with  an  extra 
unloaded  portion  of  length  £/4.  The  unloaded  portion  can  then 
be  disregarded. 

From  the  theory  of  beams,  it  is  known  that  by  applying 
external  moments  on  a beam,  the  distribution  of  the  bending 
moment  of  the  beam  can  be  changed.  It  is  therefore  suggested 
that  external  moments  be  applied  at  certain  locations  of  the 
clamped-clamped  beam  to  simulate  a double  cantilever  beam. 
However,  since  adding  external  moments  will  introduce 
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reactions  at  the  clamped  ends,  an  iterative  scheme  will  be 
needed  to  adjust  the  distribution  of  the  bending  moment. 

In  this  case,  in  order  to  make  both  of  the  loaded 
portions  of  the  clamped-clamped  beam  act  like  cantilever 
beams,  external  moments  of  the  same  amount  as  the  bending 
moments  at  x=£/4  and  x=3£/4  (calculated  when  no  external 
moment  was  added,  i.e.,  the  bending  moment  labeled  "iteration 
1"  in  figure  5.7)  are  superimposed  on  the  clamped-clamped  beam 
at  locations  just  inside  the  unloaded  portion  but  very  close 
to  x=£/4  and  x=3£/4.  When  these  moments  are  applied,  the 
resulting  bending  moments  at  x=£/4  and  x=3£/4  will  be  smaller 
than  before.  Now  this  smaller  amount  can  be  superimposed  at 
the  same  locations  again.  This  will  further  reduce  the 
resulting  bending  moments  at  x=£/4  and  x=3£/4.  By  repeating 
this  procedure,  the  bending  moment  will  finally  vanish  at 
x=£/4  and  x=3£/4  and  the  same  bending  moment  distribution  as 
that  of  a cantilever  beam  acted  upon  by  uniform  load  can  be 
obtained. 

The  process  of  iteration  is  shown  in  figures  5.6  and  5.7. 
It  can  be  seen  from  figure  5.7  that  the  moment  distribution 
in  the  loaded  region  is  gradually  reduced,  and  that  inside  the 
unloaded  region  is  gradually  increasing.  Since  the  unloaded 
region  represents  a "hole,"  the  increase  in  the  bending  moment 
inside  the  unloaded  portion  has  no  meaning.  Due  to  the 
applied  moment,  the  moment  diagram  has  a discontinuity  at 
x=i/4  and  x=3£/4. 
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Figure  5.6  Deflection  curves  of  the  double  cantilever  beam 
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The  final  result  in  both  figures  is  labeled  iteration  5. 
After  five  iterations,  the  increase  in  deflection  at  x=£/4  is 
less  than  2%  between  two  successive  iterations.  From  figure 
5.6,  it  can  be  seen  that  the  deflection  at  x=£/4  converges  to 
the  deflection  of  a uniformly  loaded  cantilever  beam  with 
length  egual  to  £/4  (the  curve  labeled  double  cantilever 
beam)  . From  figure  5.7,  it  can  be  further  understood  that  the 
moment  distribution  converges  to  that  of  a cantilever  beam. 

The  same  approach  can  now  be  applied  to  the  case  of 
plates  with  free  holes.  This  is  illustrated  in  the  following 
section  by  the  case  of  a sguare  plate  with  a central  sguare 
hole. 


Square  Plates  with  Free  Scruare  Holes 

Consider  a simply-supported  plate  which  has  a centered 
sguare  hole  as  shown  in  figure  5.8.  The  plate  is  uniformly 
loaded  with  the  load  intensity  Pg.  To  begin,  a simply- 
supported  plate  without  a hole  is  considered.  The  plate  is 
uniformly  loaded  except  for  the  portion  where  the  hole  will 
be  placed.  By  choosing  the  same  deflection  function  as 
eguation  (4.12),  the  same  formulation  for  the  strain  energy 
as  eguation  (4.13)  will  be  used. 

The  work  done  by  the  external  load  is 

W = Wp  - Wh  (5.2) 
where  Wp  is  the  work  done  by  the  external  load  when  a plate 
without  a hole  is  uniformly  loaded,  and  is  the  work  done  by 
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the  external  load  distributed  over  the  region  where  the  hole 
is  supposed  to  be  located.  In  equation  (5.2), 


W„  = 


Pq^x^v  , l-cos(mTr)^  ^ 1-cos  (n^r) 


m=l  n=l 


S S a ( 

"i  ^mn  V 


m 


) ( 


n 


) (5.3) 


and 


Wu 


Po^x^y  Nm  Nn  COS  (m7r/4 ) -COS  ( 3m7T/4  ) 


7T 


S S a r ( 

rzi  mn  L \ 


m=l  n=l 
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cos  (n7r/4)  -cos  (3n7r/4) 

( —)^  (5.4) 
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Figure  5.8  A plate  with  a free  hole 


By  using  the  work  formula  (5.2)  and  following  the  same 
procedure  as  that  shown  in  Chapter  4 , the  deflection  surface 
of  a simply-supported  plate  acted  upon  by  a uniform  load  will 
be  obtained  with  the  hole  region  not  loaded.  The  moments 
along  the  hole  edges  can  be  calculated.  Since  the  moment  of 
a plate  at  a free  edge  is  zero,  the  moments  along  the  edges 
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of  the  hole  should  be  reduced  to  zero.  To  simulate  the  effect 
of  a hole,  external  moments  are  superimposed  at  constraint 
points  chosen  from  the  edges  of  the  hole  so  as  to  reduce  the 
moment  at  each  constraint  point.  The  amount  of  external 
moment  applied  at  each  constraint  point  is  the  amount 
calculated  from  the  original  plate  (when  no  external  moment 
was  added) . As  illustrated  in  the  previous  section,  by 
iteratively  superimposing  moments  at  the  constraint  points, 
the  bending  moment  along  the  hole  edges  will  be  gradually 
reduced.  Since  the  bending  moment  in  a plate  is  continuously 
distributed  along  a section,  it  is  necessary  to  multiply  this 
moment  value  by  the  distance  between  two  successive  constraint 
points  of  each  section  in  order  to  transfer  a distributed 
moment  into  several  discrete  concentrated  moments. 

For  the  case  shown  in  figure  5.8,  the  maximum  deflection 
occurs  approximately  at  the  point  x=£^4  and  y=£„/2  or  its 
image  points  distributed  about  the  centroid  of  the  plate.  The 
maximum  deflection,  obtained  with  a 20-term  approximation,  is 
3 . 1213x10"^  (Pq£j(Vd)  for  the  sixth  iteration  and 
3 . 24080x10"^  (Pq£j^Vd)  for  the  seventh  iteration.  For  every 
iteration,  the  20-term  approximation  gives  a convergent 
deflection  surface  with  the  loading  of  that  iteration.  The 
difference  between  the  sixth  and  the  seventh  iteration  is  less 
than  4%.  The  iteration  can  be  considered  convergent. 

Figures  5.9  and  5.10  show  comparisons  between  the  finite 
element  method  with  300  elements  and  the  proposed  method  with 
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Figure  5.9  Deflection  curves  of  a plate  with  a free  hole 
taken  along  the  section  y=£  /2 
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Figure  5.10  Deflection  curves  of  a plate  with  a free  hole 
taken  along  the  section  y=£„/4 
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the  20-tenn  series  and  seven  iterations.  The  numerical  value 
of  the  maximum  deflection  calculated  from  the  finite  element 
method  with  300  elements  is  3 . 183626x10"^  (Po^^/D)  . The 
difference  between  the  finite  element  method  and  the  proposed 
method  is  less  than  2%  for  the  seventh  iteration.  The 
deflections  in  both  figures  have  been  normalized  by  the 
expression  Pq£^Vd.  The  curve  shown  in  figure  5.9  is  taken 
along  the  section  y=i^/2,  up  to  the  edge  of  the  hole.  In 
figure  5.10,  the  curve  is  taken  along  the  whole  section 
y=£y/4.  From  these  figures,  it  can  be  seen  that  the 
approximation  is  very  good. 

Figure  5.11  shows  the  deflection  surface  of  the  case 
shown  in  figure  5.8.  The  deflection  surface  has  been 
normalized  in  the  same  manner  as  before,  and  has  been 
magnified  800  times  with  respect  to  the  plate  dimensions. 


Figure  5.11 


The  deflection  surface  of  a plate 
with  a free  hole 


CHAPTER  6 

APPLICATIONS  TO  SHELL  PROBLEMS 
Introduction 

The  behavior  of  a plate  element  under  a transverse  load 
is  usually  decided  by  the  action  of  moments,  including  both 
bending  moments  and  twisting  moments,  while  a thin  shell 
generally  transmits  the  external  load  by  a membrane  stress 
which  is  uniformly  distributed  over  the  cross  section  of  the 
shell  and  normal  to  it.  This  stress  is  basically  determined 
by  the  static  equilibrium  of  the  shell  and  is  independent  of 
the  bending  effect [3].  Bending  only  has  effect  on  areas  close 
to  boundaries  and  reinforcements.  Although  the  bending  has 
only  a localized  character,  it  is  important  in  obtaining 
actual  deformations  near  the  boundaries  and  the 
reinforcements,  especially  when  the  shell  is  reinforced  with 
stiffeners  that  are  close  to  each  other.  In  view  of  the  fact 
that  the  membrane  stress  is  generally  independent  of  the 
bending  effect,  the  stresses  due  to  bending  can  be 
superimposed  on  the  membrane  stress  to  obtain  a more  accurate 
stress  distribution  in  the  region  close  to  a boundary  or  a 
reinforced  area. 
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The  bending  of  a shell  is  more  or  less  similar  to  that 
of  a plate,  and  the  theory  of  plates  may  be  applied  to  shell 
problems  in  dealing  with  bending.  Thus,  it  is  possible  to 
extend  the  approach  of  the  proposed  method  presented  in  the 
previous  chapters  to  shell  problems.  However,  it  should  be 
noted  that  the  proposed  method  can  only  be  used  to  determine 
the  bending  effect,  and  does  not  supply  any  information  about 
the  membrane  stress.  In  the  following  study,  the  focus  shall 
be  confined  to  the  linear  bending  phenomena  of  axisymmetrical 
cylindrical  shell  problems. 

For  axisymmetrical  cylindrical  shells,  it  is  possible  to 
treat  the  problem  as  a beam  supported  on  an  elastic 
foundation [ 35 ] . In  this  case,  the  modulus  of  the  "foundation" 
is  given  by 

= Eh  / r^  (6.1) 
where  E is  Young's  modulus  of  the  shell,  r is  the  radius  of 
the  cylindrical  shell,  and  h is  the  thickness.  If  the 
cylinder  is  regarded  as  a combination  of  many  axial  strips, 
then  the  effect  of  the  "elastic  foundation"  arises  from  the 
fact  that  if  any  strip  is  deflected  or  moved  radially,  the 
adjacent  strips  will  constrain  its  deflection  or  movement. 
Each  strip  will  store  an  additional  portion  of  strain  energy 
which  is  the  energy  absorbed  by  the  "elastic  foundation."  In 
other  words,  the  cylindrical  shell  can  store  more  strain 
energy  by  being  deflected  less. 
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Cylindrical  Shells  of  Infinite 
Length  Supported  bv  a Rigid  Ring 

As  mentioned  above,  bending  has  only  a localized 
character  for  shells.  The  effect  of  bending  on  the 
deflections  of  shells  reduces  rapidly  as  the  distance  from 
boundaries  or  reinforcements  increases.  This  phenomena  will 
be  illustrated  in  the  following  example.  Consider  the 
situation  shown  in  figure  6.1  where  a cylindrical  shell  is 
reinforced  or  supported  by  a rigid  circumferential  ring  (in 
order  to  have  a perfect  rigid  ring,  the  ring  is  considered 
very  thick  in  the  radial  direction) . 


V 


Figure  6.1  A rigidly  supported  cylindrical  sell 


The  coordinate  system  is  defined  as  follows  (figure  6.2)  : 
The  axial  direction  is  the  x axis,  the  radial  deflection  is 
the  V axis,  and  the  origin  is  located  at  a point  on  the  inside 
surface  of  the  shell.  With  the  coordinates  defined  as  above, 
the  deflection  v of  the  shell  is  a function  of  the  x 
coordinate  and  the  circumferential  distance  from  the  origin. 
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which  will  be  defined  as  the  circumferential  coordinate  and 
denoted  by  y. 


Figure  6.2  Coordinate  system  for  the  shell  problems 

Under  a special  circumstance,  when  both  the  geometrical 
shape  and  the  deflection  of  the  shell  are  axisymmetrical , the 
deflection  itself  is  independent  of  the  circumferential 
coordinate  y and,  thus,  can  be  considered  one-dimensional. 
Figure  6.1  represents  one  such  case.  If  the  origin  is  chosen 
on  the  rigid  support,  then  clearly  v=0  at  x=0  for  the 
situation  shown  in  figure  6.1. 

As  in  the  previous  chapters,  a deflection  function  is 
assumed  first.  Since  the  radial  deflection  of  any  point  in 
a circumferential  element  of  the  shell  is  the  same,  the 
deflection  function  is  chosen  as  follows: 

N 

v(x)  = ag  + a^cos  (n7TX/£)  (6.2) 

The  symbol  £ denotes  a reference  length.  The  reference  length 
should  be  chosen  such  that 

)S  = £[Eh/(4Dr^)  > 5 


(6.3) 
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in  order  that  the  location  at  x=£  can  be  regarded  as 
infinitely  far  away  from  the  ring,  and  bending  has  no  effect 
on  the  shell  beyond  this  location[35] . The  symbol  D was 
defined  in  eguation  (4.7).  The  strain  energy  due  to  bending 
of  the  shell  is  expressed  as 


Ui  = 


D d^v 
2 ^ dx^ 


dx  = 


2„4 


2 a/n 


n=l 


(6.4) 


J 0 

and  the  "hoop  strain  energy,"  which  is  defined  as  the  strain 
energy  stored  in  the  shell  due  to  the  effect  of  the  "elastic 
foundation,"  is 


rt- 


J 0 


2 


dx 


Eh£ 

2r^ 


+ 


N 2 

S a 

n=l  n 


(6.5) 


If  the  pressure  inside  the  cylindrical  shell  is  Pq,  then  the 
work  done  by  this  pressure  in  radially  expanding  the  shell  is 


W 


PqV(x)  dx  = Poagf 


(6.6) 


J 0 

As  before,  the  potential  energy  is 

n=Ui+U2-W  (6.7) 

Since  the  only  boundary  condition  is  that  the  radial 
defection  at  x=0  is  zero,  the  only  constraint  equation  is 


Gi(an) 


N 

a„  + 2 a„  = 0 

0 n=l 


(6.8) 
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By  following  the  same  procedure  as  described  in  the 
previous  chapters,  the  Lagrangian  function  is  formed  as 
follows: 

L = Uj  + U2  - W + (6.9) 
and  the  following  relationships  exist: 

ao  = ( Po^  - ^1  ) / ( ) (6.10.1) 

= -Aj  / n=l,2,...,N  (6.10.2) 

where 

= [ (D7rV)/(2£^)  + Y.yli/2  ] (6.11) 
By  substituting  these  equations  back  into  the  constraint 
eguation  (6.8),  the  undetermined  multiplier  Aj  can  be 
calculated.  By  substituting  the  multiplier  into  eguations 
(6.10.1)  and  (6.10.2),  the  displacement  parameters  ag  and  a^ 
(n=l, 2 , . . . ,N)  can  be  obtained. 

The  results  are  shown  in  figures  6.3  and  6.4.  Figure  6.3 
shows  the  deflection  curve  of  three  approximations  with  N=7, 
12,  and  17  respectively,  along  with  the  theoretical  deflection 
curve [35].  In  this  example,  h=0.01,  £=1.0,  pg=1.0,  r=1.0, 
E=3xl0^,  and  v=0.3.  Any  appropriate  unit  system  can  be  chosen 
for  the  above  information.  From  the  plot,  one  can  see  that 
only  ten  terms  are  needed  for  good  accuracy.  The  value  of  )3 
is  about  ten  in  this  case  and  thus,  when  x > £/2,  the  radial 
deflection  of  the  reinforced  shell  becomes  the  same  as  the 
increase  in  the  radius  of  a cylindrical  shell  which  is  free 
to  expand  without  any  constraint.  This  can  be  readily  seen 
from  the  figure.  It  can  then  be  understood  that  the  bending 
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Figure  6.3  Deflection  curves  of  a rigidly  supported  shell 
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Figure  6.4  Bending  moment  curves  of  a rigidly  supported 
shell 
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effect  has  a localized  character  and  is  only  important  in  the 
region  close  to  the  boundary  or  reinforcements.  As  a matter 
of  fact,  the  deflection  due  to  bending  is  a rapidly  damped 
oscillatory  curve  with  a wave  length  = 2irl  / )8[35].  This 
can  be  seen  from  the  figure. 

Figure  6.4  shows  the  moment  diagrams  for  the  three 
different  approximations.  As  can  be  seen  from  the  diagrams, 
the  bending  moment  dies  out  very  rapidly  after  x > £/2.  Since 
the  bending  moment  is  a quantity  resulting  from  the  second 
derivative  of  the  deflection  curve,  the  fluctuation  is  much 
more  pronounced  than  that  of  the  deflection  curve.  It  is, 
however,  possible  to  smooth  out  the  fluctuation  by  increasing 
the  number  of  terms,  as  will  be  illustrated  in  a later  case. 

Cylindrical  Shells  of  Infinite 
Length  Reinforced  by  a Narrow  Rina 

The  problem  of  a shell  reinforced  with  a narrow  circular 
ring  will  be  considered.  The  solution  can  be  obtained  by 
using  the  same  strategy  described  above.  Consider  an 
infinitely  long  cylindrical  shell  reinforced  with  a ring,  as 
shown  in  figure  6.5.  The  coordinate  system  will  be  taken  as 
that  shown  in  figure  6.2  with  the  origin  located  on  the  ring. 
Since  the  shell  is  infinitely  long,  it  may  be  taken  to  be 
symmetrical  about  the  ring,  and  thus  only  half  of  the 
arrangement  will  be  taken  into  account  in  the  following 
discussion. 
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{ 

Figure  6.5  A shell  reinforced  by  a narrow  ring 

Using  the  approximation  series  in  equation  (6.2),  the 
same  formulations  for  strain  energy  as  equations  (6.4)  and 
(6.5),  and  the  same  expression  for  work  as  equation  (6.6)  are 
obtained.  Since  the  narrow  ring  will  expand  when  the  shell 
expands  due  to  the  action  of  the  pressure  pg  inside  the 
cylindrical  shell,  strain  energy  will  be  stored.  This  strain 
energy  must  be  accounted  for  in  the  formulation.  The 

relationship  between  the  increase  of  the  radius  of  a narrow 
ring  and  the  pressure  p^,  applied  inside  the  ring  is 
linear [35] [41] : 

= p^  / K,  (6.12) 

where  v^,  is  the  increase  in  the  inside  radius  of  the  ring  and 
is  defined  as 

= EA^  / (6.13) 

In  the  above  equation,  E denotes  Young's  modulus  of  the  ring 
(which  will  be  taken  as  the  same  as  that  of  the  shell  in  this 
study)  , Aj,  is  the  cross  section  area  of  the  ring,  and  r^  is  the 
inside  radius  of  the  ring.  Since  the  thickness  of  the  shell 
is  supposed  to  be  very  thin  in  comparison  to  the  other 
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dimensions,  the  radius  of  the  ring  r,,  is  considered 
approximately  the  same  as  that  of  the  shell. 

Equations  (6.12)  and  (6.13)  are  only  valid  when  the 
thickness  of  the  ring  in  the  radial  direction  is  small  in 
comparison  with  the  radius  of  the  ring,  and  the  change  in 
radius  is  also  small  in  comparison  to  the  original  radius. 
Therefore,  in  equation  (6.13),  r^,  may  be  taken  as  the  original 
radius  of  the  ring  and  assumed  constant  over  the  process  of 
the  expansion  of  the  ring.  K,,  can  thus  be  assumed  constant 
too.  The  strain  energy  stored  in  the  ring  is 


With  the  same  assumed  deflection  function  as  in  the  previous 
section,  the  total  potential  energy  of  the  system,  including 
both  the  shell  and  the  ring,  is 


The  only  difference  between  this  case  and  the  previous 
one  is  that  now  the  deflection  is  nonzero  at  x=0.  Thus,  the 
constraint  is  that  the  defection  of  the  shell  should  be 
compatible  with  the  increase  of  the  radius  of  the  ring,  i.e.. 


U,  = K,v/  / 2 


(6.14) 


n = Uj  + U2  + u,,  - w 


(6.15) 


(6.16) 


Then  the  Lagrangian  function  takes  the  form 


L = Uj  + U2  + - W + AjGj 


(6.17) 


Following  the  same  procedure  as  before,  the  same 
relationships  as  equations  (6.10.1)  and  (6.10.2)  are  obtained. 
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By  taking  the  partial  derivative  of  L with  respect  to  v^,  an 
additional  equation  is  found: 

K^v^  - = 0 (6.18) 
By  comparing  equation  (6.18)  with  equation  (6.12),  it  can  be 
seen  that  is  the  pressure  that  acts  upon  the  ring  when  the 
inside  pressure  of  the  cylindrical  shell  is  Pq.  Thus,  Xj 
controls  the  deflection  of  the  shell  at  x=0. 

Combining  equation  (6.18)  with  the  constraint  equation 
(6.16)  and  equations  (6.10),  the  displacement  parameters  ag 
and  a^  (n=l,  2 , . . . ,N)  , and  the  increase  of  the  radius  of  the 
ring,  v^,  can  be  calculated.  Figures  6.6  and  6.7  show  the 
results  of  three  different  approximations  with  the  same 
geometry  and  material  used  in  the  previous  case,  and  Aj,=0.002. 
Figure  6.6  shows  a comparison  of  deflections  between  the 
approximations  and  the  theoretical  solution [ 35 ] . For  the 
region  far  away  from  the  reinforcement,  the  deflection 
approaches  a constant  value,  which  is  the  increase  of  the 
radius  of  a similar,  nonreinforced  cylinder  acted  upon  by  the 
same  internal  pressure.  This  is  the  same  for  both  the  present 
and  the  previous  cases,  as  would  be  expected.  Figure  6.7 
shows  the  moment  diagrams  for  the  three  different 
approximations.  It  can  be  seen  that  the  moment  decays  when 
X is  large  enough.  This  phenomena  has  been  explained  in  the 
previous  section. 
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Cylindrical  Shells  of  Finite 
Length  with  Both  Ends  Clamped 
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Figure  6.8  shows  a cylindrical  shell  of  finite  length 
fixed  at  both  longitudinal  ends.  The  coordinate  system  of 
this  problem  will  be  taken  as  the  same  system  shown  in  figure 
6.2,  with  the  origin  located  at  the  left  boundary  of  the 
cylinder  and  the  positive  x axis  running  towards  the  right. 


The  only  difference  between  the  present  case  and  the 
previous  cases  is  that  the  deflection  of  the  shell  is  zero  at 
the  both  x=0  and  x=£ . Here  £ is  no  longer  a reference  length, 
it  is  the  longitudinal  dimension  of  the  shell. 

If  the  same  deflection  function  as  equation  (6.2)  is 
assumed,  then  the  formulation  is  the  same  as  before,  except 
that  there  is  one  additional  constraint: 
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Figure  6.8  A shell  fixed  at  both  ends 


(6.19) 


Under  this  condition,  the  Lagrangian  function  is 


L — + U2  ~ W + + A2G2 


(6.20) 
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The  displacement  parameters  are  as  follows,  using  the  same 
procedure  illustrated  in  the  previous  sections: 

ao  = ( Po^  - - ^2  ) / ( V ) (6.21.1) 

a„  = -[  Aj  + A2COs(n7T)  ] / (6.21.2) 
Substituting  these  relationships  (6.21.1)  and  (6.21.2)  into 
the  constraint  equations  (6.8)  and  (6.19)  will  result  in  two 
linear  equations  with  the  undetermined  multipliers  Aj  and  A2 
as  unknowns.  It  then  follows  that  the  displacement  parameters 
can  be  obtained  by  using  equation  (6.21). 

Using  the  same  geometry  and  material  as  in  the  previous 
cases,  the  results  shown  in  figures  6.9  and  6.10  are  obtained. 
In  figure  6.9,  the  deflections  of  the  shell  with  three 
different  approximations,  N=7,  12,  and  17  respectively,  are 
shown.  The  theoretical  deflection  derived  from  the  basic 
theory  of  shells  is  also  shown  for  comparison.  Since  the 
difference  between  the  12-term  approximation  and  the 
theoretical  curve  is  almost  neligible,  the  12-term  series  may 
be  taken  as  a very  good  approximation  of  the  deflection. 
Again,  since  /3  (referring  to  equation  6.3)  is  about  twice 
five,  the  deflection  of  the  shell  at  the  region  around  x=£/2 
is  approximately  the  same  as  a similar  cylindrical  shell 
without  reinforcements.  Figure  6.10  shows  the  moment  diagrams 
corresponding  to  the  approximate  deflection  curves  shown  in 
figure  6.9.  The  moment  is  approximately  zero  around  x=£/2  for 


the  same  reason. 
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Figure  6.9  Deflection  curves  of  a shell  with  both  ends 
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Figure  6.10  Bending  moment  curves  of  a shell  with  both 
ends  clamped 


Clamped  Cylindrical  Shells  with 
Circumferential  Ring  Reinforcements 
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The  cases  discussed  so  far  can  be  further  extended  to  a 
cylindrical  shell  clamped  at  both  ends  and  reinforced  with 
rings  as  shown  in  figure  6.11.  The  coordinate  system  will 
remain  the  same. 


Figure  6.11  A clamped  shell  reinforced  by  rings 

The  strain  energy  stored  in  each  ring  should  be  added  to 
the  total  potential  energy,  i.e., 

n = Uj  + U2  + S - W (6.22) 
where  is  defined  the  same  as  in  equation  (6.14),  with  the 
subscript  i denoting  the  ith  ring.  The  summation  in  the  above 
equation  should  be  taken  over  all  rings. 

The  compatibility  constraint  states  that  the  deflection 
of  the  shell  should  be  compatible  with  the  radial  expansion 
of  the  rings.  Besides  equations  (6.8)  and  (6.19),  the 
following  equations,  which  define  this  compatibility,  should 
also  be  part  of  the  constraint  equations: 
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G3i(an)  = ( ao  + a,cos(n7Tx/£)  ) - = 0 

(6.23) 

where  is  the  radial  expansion  of  the  ith  ring  and  x^  is  its 
location  in  the  x axis.  is  the  number  of  rings.  By 

following  the  same  procedure  outlined  in  the  previous 
sections,  the  following  equations  for  the  displacement 
parameters  can  be  obtained: 

ao  = ( Po^  - - ^2  - |i/^i  ) / ( ) (6.24.1) 

an  = -[  ^1  + A2Cos(n7T)  + 

i=l,2,...,I^  (6.24.2) 

In  equations  (6.24.1)  and  (6.24.2),  and  A2  are  the 

undetermined  multipliers  associated  with  the  clamped  boundary 
conditions  (6.8)  and  (6.19),  and  (i=l,  2 , . . . , I,,)  are 

associated  with  the  deflection  compatibility  required  by 
equation  (6.23).  Also,  the  following  equations  for  the  radial 
expansion  of  the  rings  are  found: 

= Mi  / K^i  i=l,2,...,I^  (6.24.3) 

Substituting  these  equations  into  the  constraint 
equations  (6.8),  (6.19),  and  (6.23)  results  in  a set  of  (1^+2) 
simultaneous  linear  equations  with  the  undetermined 
multipliers  A^,  A2,  and  Mi  (i=l/ 2 , . . . , I,,)  as  unknowns.  Solving 
these  simultaneous  equations  and  substituting  the  solution 
back  into  equation  (6.24)  gives  the  displacement  parameters. 
The  radial  deflection  of  the  shell  can  then  be  calculated. 
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With  equal  to  two  and  the  rings  equally  spaced  on  the 
shell  (figure  6.11),  the  result  is  shown  in  figure  6.12.  The 
geometry  and  material  of  the  shell  are  the  same  as  in  the 
previous  cases.  For  the  deflection  curve  shown  in  figure 
6.12,  a 17-term  approximation  was  sufficient  to  obtain 
convergence.  Since  it  is  extremely  difficult  to  fomnulate  the 
theoretical  deflection  for  this  case,  it  is  not  available  and 
thus  not  shown  in  the  figure  for  comparison.  Figure  6.13 
shows  the  moment  diagram  which  is  obtained  from  a 100-term 
approximation.  The  diagram  is  much  smoother  than  the  moment 
diagrams  for  the  previous  cases.  This  is  because  more  terms 
are  taken  to  approximate  the  actual  deflection,  and  thus  the 
deflection  curve  is  in  excellent  agreement  with  the  actual 
deflection  curve.  It  should  be  noted  that  100  terms  is  much 
more  than  necessary.  Thirty  terms  will  give  smooth  curves, 
even  for  the  moment  diagram. 

Shells  Reinforced  by  Rings 
and  Longitudinal  Stiffeners 

As  a general  application  to  shell  problem,  the  method  of 
treating  a cylindrical  shell  reinforced  by  both 
circumferential  rings  and  longitudinal  stiffeners  will  be 
presented  in  this  section.  Figure  6.14  shows  a general 
arrangement  of  how  shells  can  be  reinforced  by  both  rings  and 
longitudinal  stiffeners.  In  this  situation,  the  deflection 
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of  the  shell  is  no  longer  a one-dimensional  problem,  and  thus 
no  theoretical  solution  is  available. 


Figure  6.14  A shell  reinforced  by  rings  and 
longitudinal  stiffeners 

Using  the  same  coordinate  system  as  before,  the 
deflection  function  of  the  shell  can  be  taken  as 

Nm  Nn 

v(x,y)  = S S a„cos(2m7rx/£)cos(ny)  (6.25) 

m=0  n=0 

The  deflection  functions  of  the  longitudinal  stiffeners  are 
given  by 

Nm 

^bi(x)  = ^oCn,i(  1-COS  (2m7rx/£)  ) 

i=l , 2 , . . . , (6.26) 

where  1,^  is  the  number  of  the  longitudinal  stiffeners.  The 
deflection  functions  of  the  rings  (expansions  of  the  rings) 
are 

Nn 

Vrj(x)  = S^d^jCos(ny)  j=l,2,...,I^  (6.27) 

where  is  the  number  of  rings.  If  the  placement  of  the 
reinforcements  is  symmetric  about  the  center  of  the  cylinder. 
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then  only  a quarter  of  the  shell  needs  to  be  considered.  In 
this  case  the  bending  strain  energy  of  the  shell  is  given  by 


ptn  pj 


Ui  = (Dr/2) 


0 


( + 3^v/5y^  dx  dy 


0 


And  the  hoop  strain  energy  of  the  shell  is 

ptn  pf 


U2  = (Khr/2) 


yr  dx  dy 


(6.28) 


(6.29) 


The  total  potential  energy  stored  in  the  shell  is 


U = + U2 


Nm  Nn  9 

S S a a 

m=0  n=0 


(6.30) 


where 


“mn  = Kh7rr£/4  if  m=n=0 

= (K^,7rr£/8)  + (D7rr£/8 ) (2m7r/£ ) "^ 

if  iii^O  and  n=0 
= (Kp,7rr£/8)  + (DTrrfnVs) 


if  m=0  and  n^O 

= (Kh7rr£/16)  + (D7rr£/16)  [ (2iti7r/£)  (6.31) 

if  m#0  and  n^O 

The  bending  energy  stored  in  each  longitudinal  stiffener  is 

ptn 


u 


bi 


(EI/2) 


(3\/3x^)^  dx 


J 0 


i=l,2, 


(6.32) 


where 
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Tn,  = (EI«/8)  (2m7r/£) 


(6.33) 


The  strain  energy  stored  in  each  ring  is 


= (K,r/2) 


Vrj  dy 


Nn  9 
= S d . 6 
n=0  "J  " 


1 1,2, ...,X 


(6.34) 


where 

5,  = K,r7T/2 
= K^rTT/4 

The  work  done  by  the  pressure  Pg  is 

n«/2  rf 


if  n=0 
otherwise 


W = r 


PoV(x,y)  dx  dy 


(6.35) 


Nm  Nn 

= S 2 a„„n„„ 

m=0  n=0 

where 


(6.36) 


= ^l=-«Po/2 


if  m=n=0 


(6.37) 


= 0 otherwise 

The  constraints  are  the  zero  deflection  at  the  end  x=0 
and  the  deflection  compatibility  between  the  shell,  the 
stiffeners,  and  the  rings.  These  are 

v(0,  y,^)  = 0 k=l,2,  . . . ,Kbi  (6.38.1) 


v(x^,  y.)  - Vj,^(Xk)  = 0 

i=l,2, ,1^,,  k=l,2, (6.38.2) 

where  denotes  the  number  of  constraint  points  along  the 
boundary  x=0,  and  denotes  the  number  of  constraint  points 
chosen  from  the  ith  longitudinal  stiffener.  Also, 
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v(Xj,  Yk)  - v,j(Yk)  = 0 

j=l,  2 , . . . , I^,  k=l,2,  . . . (6.38.3) 

where  K^j  denotes  the  number  of  constraint  points  chosen  from 
the  jth  ring.  In  order  to  simplifY  the  formulation,  assume 
that  the  number  of  constraint  points  on  the  boundarY  x=0  and 
all  reinforcements,  including  the  rings  and  the  stiffeners, 
are  the  same  (K^,)  , i.e.,  K(j^=K^^=Kjj=K(,.  The  constraint 

equations  can  then  be  expressed  in  the  following  compact  form: 

N N N 


s 

tti=0 


s 

n=0 


‘mm  j 


- S 

m=0 


= 0 


j=l,2,...,K^  (6.39) 


If  the  constraint  points  are  regarded  as  points  taken  from 
"constraint  lines"  which  are  on  the  surface  of  the  CYlindrical 
shell,  then  the  subscript  i denotes  the  ith  constraint  line. 
These  constraint  lines  include  the  boundarY  x=0,  the 
longitudinal  stiffeners,  and  the  rings.  The  subscript  j 
denotes  the  jth  constraint  point  on  the  ith  constraint  line. 
The  definition  of  R„,j,  and  f^^  are  listed  as  follows: 

1.  When  the  constraint  point  is  on  the  boundarY  ^=0 

Qmnij  = cos(nYj) 

Rmj  = 0 
fmi  = 0 

2 . When  the  constraint  point  is  on  a longitudinal 
stiffener 


Qmnij  = cos(2m7rXj/€)cos(nYi) 
= 1 - cos(2m7rx/£) 

^mi  ~ *^mi 
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3 . When  the  constraint  point  is  on  a ring 
Qmnij  = cos(2m7rx/£)cos(nyj) 

Rn,j  = cos(nYj) 
fn,i  = 

The  following  eguations  relating  the  displacement  parameters 
a,^^  and  to  the  undetermined  multipliers  can  then  be 
derived: 


Ic  Kc 


(20=n,n) 

(6.40) 

and 

f.l  = 

Kc 

J=1 

) / (2CJ 

(6.41) 

where 

Tn, 

when  f^,  = c^, 

(6.42.1) 

= 

<5. 

when  f^^  = d^^ 

(6.42.2) 

and 

Ic  = 

Ib  + 

Ir  + 1 

(6.43) 

In  the 

above 

equation,  I^,  denotes 

the  total 

number  of 

constraint  lines.  By  substituting  equations  (6.40)  and  (6.41) 
into  equation  (6.39),  a set  of  linear  equations  can  be 
obtained.  Solving  these  equations  gives  the  undetermined 
multipliers  A.j.  Then,  by  substituting  the  undetermined 
multipliers  into  equations  (6.40)  and  (6.41),  the  parameters 
a^j^  can  be  calculated,  and  the  deflection  surface  can  be  found. 

Figure  6.15  shows  the  circumferential  deflection  surface 
of  a shell  with  length  0.7  units,  radius  1.0  unit,  and 
thickness  0.01  units.  The  cylindrical  vessel  is  reinforced 
by  four  rings,  and  four  longitudinal  stiffeners,  with 
thickness  0.02  units  and  cross  section  height  equal  to  0.04 
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units.  The  result  is  obtained  by  assuming  a 25-term 
approximation  function  (N^=N^=25)  as  given  in  eguation  (6.25). 
This  number  of  terms  is  large  enough  to  provide  a convergent 
result.  The  trace  of  the  reinforcements  can  be  identified  in 
the  drawing.  The  deflection  of  the  circumferential  surface 
drawing  has  been  magnified  27,000  times  relative  to  the 
dimensions  of  the  shell  in  order  to  provide  detail. 


O 


Figure  6.15  The  circumferential  deflection  surface  of 
a shell  reinforced  by  rings  and 
longitudinal  stiffeners 


CHAPTER  7 

APPLICATIONS  TO  STRUCTURAL  DESIGN 

The  previous  chapters  illustrate  the  use  of  the  Rayleigh- 
Ritz  method  incorporating  undetermined  multipliers  to  the 
analysis  of  structural  elements  such  as  beams,  plates,  and 
shells.  As  was  seen,  the  proposed  method  provides  an  elegant, 
efficient,  and  effective  way  for  the  analysis  of  elastic 
structural  elements. 

The  proposed  method  can  also  be  used  as  a valuable  design 
tool,  especially  for  complicated  structures,  such  as  the 
honeycomb  discussed  in  Chapter  4 . Designers  so  far  almost 
exclusively  rely  on  the  finite  element  method  to  help  in  the 
design  of  this  type  of  complicated  structure.  In  view  of  the 
great  effort  of  preprocessing  and  postprocessing  that  the 
finite  element  method  requires,  however,  the  method  can  be 
very  expensive  considering  the  iterative  nature  of  optimum 
design.  On  the  other  hand,  by  taking  advantage  of  the 
modified  Rayleigh-Ritz  method,  the  optimum  design  of  complex 
structural  elements  can  be  carried  out  with  relative  ease. 
The  proposed  method  can  save  a great  amount  of  time  and  labor 
used  in  the  iterations  required  in  the  optimum  design  process. 
In  the  following  sections,  some  optimum  design  cases  are 
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illustrated  using  the  modified  Rayleigh-Ritz  method  for  the 
analytical  formulation. 

Beam  Design 

Since  the  proposed  method  provides  a convenient  way  to 
analyze  the  problem  of  multi-span  beams,  it  is  can  be  used  as 
a design  tool  for  this  type  of  structure.  Consider  the  case 
that  was  shown  in  figure  3.12,  with  the  span  length,  or  the 
location  of  the  interior  support  x^,  assumed  to  be  variable. 
Since  is  variable,  it  is  possible  to  make  the  deflection 
curve  of  the  beam  as  close  as  possible  to  some  prespecified 
form  by  adjusting  x^.  An  "error"  function,  which  indicates 
the  "closeness"  of  the  approximate  deflection  curve  obtained 
from  the  modified  Rayleigh-Ritz  method  to  the  prespecified 
points  can  therefore  be  defined  in  the  following  form: 

R = ( v^(Xj)  - Uj  (7.1) 

where  R denotes  the  "error,"  Uj  denotes  the  prespecified 
deflections  at  the  particular  locations  x=Xj  ( j=l, 2 , . . . , J) , 
and  J is  the  number  of  prespecified  points.  The  approximation 
deflection  function  v^^  is  taken  to  be  the  same  as  equation 
(3.22). 

* 

The  next  step  is  to  find  the  particular  x^  which 
minimizes  R.  (Remember,  for  this  case,  R is  a function  of  the 
location  of  the  interior  support  only.)  Since  there  is  only 
one  independent  variable  x^,  the  Newton-Raphson  method [42]  is 
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an  effective  method  to  find  the  minimum  point  . It  should 
be  noted  that  for  any  specific  x^,  the  displacement  parameters 
ag,  a„,  and  (n=l, 2 , . . . ,N)  of  the  assumed  deflection  function 
v^(x)  can  be  calculated  with  the  modified  Rayleigh-Ritz 
method,  and  then  the  error  R for  that  specific  x^  can  be 
computed.  The  first  and  the  second  derivatives  of  R with 
respect  to  x^  are  approximately  calculated  as  follows: 

dR/dx^  = ( RCx^+Ax^)  - RCXg)  ) / AX3  (7.2) 

and 

d^R/dx3^  = ( R(X3+2ax3)  - 2R(X3+ax3)  + R(X3)  ) / (ax3)^ 

(7.3) 

where  AXg  is  a small  increment  in  X3.  Once  again,  the 
displacement  parameters  for  the  locations  X3+AX3  and  X3+2AX3 
are  obtained  with  the  modified  Rayleigh-Ritz  method.  With 
these  parameters,  the  errors  R(x3+ax3)  and  R(X3+2ax3)  can  then 
be  computed. 

Using  the  prespecified  deflections  outlined  in  table  7.1, 
the  Newton-Raphson  method,  incorporating  the  modified 
Rayleigh-Ritz  method,  gives  the  result  shown  in  figure  7.1  in 
five  iterations. 


Table  7.1  Specification  of  deflection 


0.1 

-0.00035 

0.3 

0.00055 

0.5 

0.00330 

5.0e  — 003 
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In  the  figure,  the  symbol  X represents  the  desired 
deflection.  The  highest  curve  corresponds  to  the  result  of 
the  first  iteration  and  the  other  two  represent  the  third  and 
fifth  iterations  using  the  Newton-Raphson  scheme, 
respectively.  The  curve  of  the  fifth  iteration  (the  final 
configuration)  is  almost  coincident  with  all  the  X's.  Note 
that,  in  general,  the  final  configuration  is  not  always  close 
to  the  specified  curve,  and  how  close  the  approximate 
deflection  curve  will  be  to  the  prespecified  deflection  is 
dependent  upon  the  shape  of  the  specified  curve. 

Similarly,  an  optimum  location  for  an  external 
concentrated  force  can  be  found  such  that  the  deflection  curve 
will  be  as  close  as  possible  to  some  prespecified  points. 
Taking  a simply-supported  beam  as  an  example  (figure  3.1)  with 
the  distributed  load  p(x)  replaced  by  a concentrated  load,  and 
specifying  the  required  deflection  condition  (Uj)  as  outlined 
in  table  7.2,  the  result  obtained  from  seven  iterations  of  the 
Newton-Raphson  method,  using  the  modified  Rayleigh-Ritz 
procedure  in  each  iteration  to  calculate  the  displacement 
parameters,  is  shown  in  figure  7.2.  Symbol  X denotes  the 
desired  deflections  at  the  locations  Xj  ( j=l, 2 , . . . , J) . It 
should  be  understood  that  the  general  procedure  is  similar  to 
the  previous  one,  except  that  the  independent  variable  is  now 
the  location  of  the  load,  not  that  of  the  interior  support. 

Since,  in  this  study,  the  deformation  of  a body  is 
confined  to  be  within  the  linearly  elastic  range,  the 
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deformations  that  a body  undergoes  due  to  several  loads  can 
be  superimposed.  If  a simply-supported  beam  is  acted  upon  by 
more  than  one  concentrated  load,  such  as  the  one  shown  in 
figure  7.3  where  two  concentrated  loads  are  applied  to  the 
beam,  the  final  configuration  that  the  beam  takes  is  simply 
the  sum  of  the  configurations  due  to  each  individual  load. 


where  denotes  the  number  of  concentrated  loads,  denotes 
the  magnitude  of  the  ith  load,  and  u^(x)  is  the  deflection 
curve  due  to  a unit  concentrated  load  applied  at  the  same 
location  as  the  ith  load.  These  deflection  curves  due  to  unit 
loads  can  be  obtained  by  the  modified  Rayleigh-Ritz  method. 


i.e.  , 


(7.4) 


Table  7.2  Specification  of  defection 


0.2 

0.4 

0.6 

0.8 


0.01200 

0.01920 

0.01813 

0.01067 


In  this  case,  the  error  function  is  defined  as 


(7.5) 


In  the  least  square  sense,  equation  (7.5)  represents  the 
criterion  for  fitting  the  discrete  data  point  Uj  ( j=l, 2 , . . . , J) 
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with  the  functions  ( i=l , 2 , . . . , ) . It  is  relatively  easy 

to  obtain  the  magnitude  for  each  load  (assuming  that  the 
location  of  is  known)  so  that  the  error  is  minimum.  This 
can  be  done  by  taking  partial  derivatives  of  R with  respect 
to  each  P^,  setting  each  derivative  equal  to  zero,  and 
simultaneously  solving  the  equations.  The  equations  take  the 
following  form: 

jli  ( li 

k=l,2,...,If  (V.6) 


■XiH 


V 


X2- 


X 


Figure  7.3  A simply-supported  beam  acted  upon  by 
two  concentrated  loads 

By  applying  these  techniques,  it  is  possible  to  simulate 
the  deflection  curve  of  a beam  due  to  an  arbitrarily 
distributed  load.  The  deflection  curve  due  to  the  arbitrary 
load  can  be  regarded  as  the  prespecified  curve.  A sinusoidal 
load  distribution  is  taken  as  an  example,  and  five 
concentrated  loads  are  used  to  simulate  the  distributed  load. 
The  locations  of  the  concentrated  loads  are  x=f/10,  3£/10, 
5£/10,  7f/10,  and  9£/10.  The  deflection  curves  due  to  unit 
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concentrated  loads  in  these  locations  can  be  obtained  using 
the  modified  Rayleigh-Ritz  method.  As  stated  in  equation 
(7.4),  these  curves,  when  respectively  multiplied  by  the 
magnitudes  P.  and  added  together,  give  the  final 
configuration.  The  resulting  deflection  v(Xj) , along  with  the 
prespecified  deflection  Uj,  is  shown  in  table  7.3.  Due  to 
the  symmetry  about  the  center  of  the  beam,  only  half  of  the 
data  is  shown  in  the  table.  Note  that  if  the  number  of 
prespecified  points  is  the  same  as  that  of  the  concentrated 
loads,  then  the  deflection  curve  can  be  forced  to  pass  through 
all  the  specified  points. 


Table  7 . 3 Result  and  comparison  of  the  case  where  the 
magnitudes  of  external  loads  are  unknown 


v(xp 

error (%) 

0.05 

1.606x10'^ 

1.607x10'^ 

0.062 

0.10 

3.172x10'^ 

3.173X10'^ 

0.032 

0.15 

4.661x10'^ 

4.660x10'^ 

0.021 

0.20 

6.034x10'^ 

6.033x10’^ 

0.017 

0.25 

7.259x10'^ 

7.260x10’^ 

0.014 

0.30 

8.305x10'^ 

8.307x10'^ 

0.024 

0.35 

9.147x10'^ 

9.147x10’^ 

0.000 

0.40 

9.763x10'^ 

9.761X10'^ 

0.020 

0.45 

1.014x10'^ 

1.014x10'^ 

0.000 

0.50 

1.027x10’^ 

1.027x10'^ 

0.000 

By  combining  the  approaches  used  above  to  determine  the 
optimum  location  and  the  optimum  magnitude  of  a concentrated 
load,  it  is  possible  to  obtain  an  optimum  load  whose  location 
and  magnitude  are  not  known.  Since  the  equation  (7.6)  used 
to  solve  for  the  magnitude  is  linear,  it  is  always  possible 


145 


to  obtain  the  best  magnitude  corresponding  to  a given  location 
directly,  without  applying  Newton-Raphson  method  which 
reguires  iterations.  This  suggests  that  it  may  be 

advantageous  to  treat  the  problem  in  an  alternate  fashion. 
In  other  words,  to  find  the  magnitude  from  the  location  and 
then  find  the  location  from  the  magnitude.  A conceptual 
algorithm  for  accomplishing  this  can  be  stated  as  follows: 

1.  Start  with  an  initial  guess  of  the  location  of  the 
load. 

2 . Calculate  the  best  magnitude  corresponding  to  the 
current  location  with  equation  (7.6). 

3.  Calculate  the  "error"  with  the  current  location  and 
magnitude . 

4.  With  the  magnitude  held  constant  temporarily,  improve 
the  location  with  the  Newton-Raphson  scheme. 

5.  Check  if  it  is  convergent — if  not,  repeat  the 
procedure. 

In  table  7.4,  the  result  of  an  example  is  shown  with  the 
iteration  process  outlined  in  table  7.5.  Each  of  the 
intermediate  results  was  obtained  by  the  modified  Rayleigh- 
Ritz  method.  The  prespecified  deflections  Uj  are  calculated 
from  the  case  of  a simply-supported  beam  acted  upon  by  a 
concentrated  load  with  P=-0.5  and  Xp=0.8£.  It  can  be  seen 
from  table  7.5  that  the  error  gradually  decreases  as  the 
iteration  proceeds.  As  shown  in  table  7.4,  the  final  result 
v(Xj)  is  almost  exactly  the  same  as  the  prespecified 
deflection.  Nevertheless,  how  small  the  minimum  error  will 
be,  or  how  close  the  final  configuration  will  be  to  the 
prespecified  points,  is  dependent  upon  the  shape  of  the 
specified  curve. 
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Table  7 . 4 Result  of  the  case  where  magnitudes 
of  external  loads  are  unknown 


Xj  Uj  V(Xj) 


0.05 

-7.979x10'" 

-7.979x10'" 

0.10 

-1.583x10'^ 

-1.583x10'^ 

0.15 

-2.344x10'^ 

-2.344x10'^ 

0.20 

-3.067x10'^ 

-3.067x10'^ 

0.25 

-3.740x10'^ 

-3.740x10'^ 

0.30 

-4.350x10'^ 

-4.350x10'^ 

0.35 

-4.885x10'^ 

-4.885X10'^ 

0.40 

-5.333x10'^ 

-5.333x10'^ 

0.45 

-5.681x10'^ 

-5.681x10'^ 

0.50 

-5.917x10'^ 

-5.917x10'^ 

0.55 

-6.027x10'^ 

-6.027x10'^ 

0.60 

-6.000x10'^ 

-6.000x10'^ 

0.65 

-5.823x10'^ 

-5.823x10'^ 

0.70 

-5.483x10'^ 

-5.483x10'^ 

0.75 

-4.969x10'^ 

-4.969x10'^ 

0.80 

-4.267x10'^ 

-4.267x10'^ 

0.85 

-3.375x10'^ 

-3 . 375x10'^ 

0.90 

-2.333x10'^ 

-2.333x10'^ 

0.95 

-1.192x10'^ 

-1.192x10'^ 

Table  7.5  Intermediate  results  of  iteration 


iteration 

location 

magnitude 

R 

1 

0.500000 

-0.293777 

2.2x10'^ 

2 

0.822525 

-0.555668 

7.6x10'^° 

3 

0.797026 

-0.494130 

1.4x10'^^ 

4 

0.799914 

-0.500245 

1.2X10'^" 

5 

0.799950 

-0.500322 

4.0x10'^^ 

The  same  approach  can  be  applied  to  solve  the  optimum 
design  problem  of  a beam  acted  upon  by  several  loads  with 
either  the  locations,  the  magnitudes,  or  both  not  specified. 
Furthermore,  by  using  the  Rayleigh-Ritz  method  together  with 
the  Lagrangian  rule  (the  modified  Rayleigh-Ritz  method) , the 
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design  of  a beam  with  any  boundary  conditions  acted  upon  by 
several  loads,  for  which  the  magnitudes  and  locations  are  not 
known,  can  be  accomplished. 

Consider  the  case  shown  in  figure  7.3  where  two 
concentrated  loads  are  applied  to  a simply-supported  beam. 
If  the  magnitudes  of  the  loads  are  known,  then  there  are  two 
decision  variables,  namely  the  two  locations  Xj  and  Xg  where 
the  loads  are  applied.  Taking  the  specifications  outlined  in 
table  7.6  as  an  example,  the  optimum  locations  of  x^=0.5£  and 
X2=0.8£,  which  result  in  an  error  of  the  order  of  10'^\  are 
obtained  within  an  average  of  35  iterations  from  different 
starting  points. 


Table  7.6  Specification  of  the  case  where  the 
locations  of  loads  are  unknown 


^3 

«3 

^3 

^3 

0.05 

2.530x10'^ 

0.55 

1.821x10'^ 

0.10 

5.017x10'^ 

0.60 

1.790x10'^ 

0.15 

7.415x10'^ 

0.65 

1.714x10'^ 

0.20 

9.683x10'^ 

0.70 

1.592x10'^ 

0.25 

1.178x10’^ 

0.75 

1.432x10'^ 

0.30 

1.365x10’^ 

0.80 

1.208x10'^ 

0.35 

1.526x10'^ 

0.85 

9.478x10’^ 

0.40 

1.657x10'^ 

0.90 

6.517x10'^ 

0.45 

1.752x10'^ 

0.95 

3.318x10'^ 

The  solution  method  is  briefly  stated  as  follows:  One 
of  the  locations,  say  the  second  one  (X2)  is  first  assumed  to 
be  given  and  fixed.  Then  the  Newton-Raphson  method  is  used 
to  obtain  a better  location  for  the  first  one,  with  the 


148 


intermediate  result  obtained  by  the  modified  Rayleigh-Ritz 
method.  Next,  the  location  of  the  first  one  is  assumed  fixed, 
and  the  Newton-Raphson  method,  incorporating  the  modified 
Rayleigh-Ritz  method,  is  applied  to  find  a better  location  for 
the  second  one.  The  process  is  repeated  until  the  optimum 
locations  are  found. 

In  comparison  with  the  procedure  stated  above,  the 
classical  optimization  schemes  have  the  disadvantages  of 
either  taking  more  iterations  to  obtain  the  solution  or 
needing  a good  guess  of  the  starting  point  in  order  to  avoid 
divergence.  An  example  of  the  former  is  the  steepest  descent 
method  [32],  and  that  of  the  later  is  the  Hooke  and  Jeeves 
method [32] . 


Plate  Design 

In  this  section,  the  application  of  the  modified 
Rayleigh-Ritz  method  to  the  analytical  aspect  of  the  optimum 
design  of  plate  structures  will  be  illustrated.  As  an 
example,  consider  the  structure  shown  in  figure  4.20,  which 
is  a honeycomb  sandwich  structure.  This  type  of  structure  can 
hardly  be  solved  by  conventional  methods.  It  is  also 
difficult  for  the  finite  element  method,  which  is  one  of  the 
most  flexible  numerical  schemes,  although  expensive  in  the 
sense  of  computer  time.  The  practicability  and  efficiency  of 
the  modified  Rayleigh-Ritz  method  makes  it  a very  good  choice 
for  the  optimum  design  of  structures  of  this  kind. 
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The  strength  and  weight  (or  equivalently,  the  volume)  of 
a structure  are  two  important  factors  which  a designer  must 
take  into  account.  For  designing  the  considered  structure, 
the  decision  variables  may  include  the  geometrical  dimensions 
of  the  overall  structure  or  those  of  the  individual  plate 
elements  (the  top  and  the  bottom)  and  the  stiffeners.  To 
simplify  the  problem,  assume  first  that  both  the  top  and 
bottom  plates  are  identical,  and  the  stiffeners,  although  they 
may  have  different  lengths,  have  identical  vertical  cross 
sections.  Further,  assume  that  the  dimensions  of  the 
structure  are  given.  The  only  decision  variables  left  are  the 
thicknesses  of  the  plate  elements  and  the  stiffeners,  provided 
the  number  and  the  locations  of  the  stiffeners  are  also  given. 

Since  the  stiffeners  are  usually  slender,  it  is 
reasonable  to  determine  the  thickness  of  the  stiffeners  in 
such  a way  that  buckling  in  the  vertical  direction  can  be 
avoided.  The  critical  buckling  load  for  a beam  can  be 
approximately  determined  from  the  following  equation[43] [44] ; 

Per  = / h3^  (7.7) 
where  h^  is  the  height  of  the  stiffeners,  which  is  equal  to 
the  overall  height  of  the  structure  (H)  minus  the  thickness 
of  both  of  the  top  and  bottom  plates  (h) , i.e.,  h^  = H - 2h. 
Ijj  is  the  area  moment  of  inertia  per  unit  length  of  the 
stiffeners,  taken  horizontally: 

Ib  = t3^  / 12 

where  t3  is  the  thickness  of  the  stiffeners. 


(7.8) 
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Assuming  that  the  yield  stress  of  the  structure  is  a^, 
then  the  maximum  uniform  load  that  a "solid  plate"  can  take 
without  exceeding  its  yield  point  is [3] 

p^  = / (0.7134£/)  (7.9) 

and  the  average  load  per  unit  length  of  the  stiffeners  is 

p=p££/£  (7.10) 

where  denotes  the  total  length  of  all  the  stiffeners.  A 
reference  "solid  plate"  is  defined  in  this  investigation  as 
a regular  plate  with  the  same  dimensions  (length,  width  and 
height)  and  same  material  as  that  of  the  honeycomb  structure. 
If  denotes  the  safety  factor  for  buckling,  then  the 

following  equation  can  be  obtained: 

Per  - Ke.f.tvP.s  (’-PI* 

The  forces  p^  and  p^^  can  be  determined  from  equations 
(7.9)  and  (7.10)  with  the  geometry  information  of  the 
structure  and  the  strength  of  the  material  given.  By 
substituting  equations  (7.7),  (7.8),  (7.9),  and  (7.10)  into 

equation  (7.11),  a relationship  between  tg  and  hg  can  be 
obtained.  When  H is  given,  hg  can  be  determined  by  h.  The 
thickness  of  the  stiffeners  for  any  particular  h can  then  be 
found  to  insure  that  buckling  is  avoided. 

Figure  7 . 4 shows  the  relationship  between  the  thicknesses 
of  the  plate  elements  and  the  stiffeners,  with  the  same 
geometry  and  material  of  the  honeycomb  structure  as  was  shown 
in  Chapter  4.  It  is  seen  from  the  figure  that  as  the 
thickness  of  the  plate  elements  increases  and  approaches  0.05, 
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which  is  half  of  the  overall  height  of  the  honeycomb 
structure,  the  thickness  of  the  stiffeners  decreases  and 
approaches  zero.  Since  the  overall  height  of  the  structure 
is  assumed  to  be  given,  it  can  be  seen  that  as  the  thickness 
of  the  plate  elements  increases,  the  height  of  the  stiffeners 
is  decreased  and,  consequently,  the  stiffeners  can  still 
support  the  same  buckling  load  with  the  smaller  thickness. 
For  any  given  thickness  of  the  plate  elements,  the  thickness 
of  the  stiffeners  can  be  determined  from  equations  (7.7  - 
7.11)  and  the  total  volume  of  the  structure  can  be  determined 
by  adding  the  volume  of  the  stiffeners  and  the  plate  elements 
together.  Figure  7.5  shows  the  change  in  total  volume  as  the 
thickness  of  the  plate  elements  increases. 

Having  determined  the  thickness  of  the  stiffeners  in  the 
manner  described  above,  the  remaining  decision  variable  is  the 
thickness  of  the  plate  elements  of  the  honeycomb  structure. 
Usually,  in  designing  a structure,  two  quantities,  strength 
and  rigidity,  can  be  chosen  as  reference.  If  strength  is  more 
important  than  rigidity,  the  issue  becomes  what  is  the  best 
thickness  of  the  plate  elements  that  gives  the  most  gain  in 
strength.  On  the  other  hand,  putting  emphasis  on  the  rigidity 
means  finding  the  best  thickness  of  the  plate  elements  that 
gives  the  most  gain  in  rigidity. 

To  have  a comparison  standard  for  strength,  the  overall 
deflection  of  the  structure  is  taken  as  the  reference  to 
construct  a plate  with  uniform  thickness  for  comparison,  which 
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Figure  7.5  Volume  of  the  honeycomb  structure 
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will  be  called  the  equivalent  plate  (whose  length  and  width 
are  the  same  as  those  of  the  honeycomb  structure)  . The 
thickness  of  the  equivalent  plate  is  so  determined  that  the 
maximum  deflection  of  the  plate  is  the  same  as  the  maximum 
overall  deflection  of  the  honeycomb  structure,  when  they  are 
acted  upon  by  the  same  load.  The  ratio  of  the  maximum  stress 
of  the  honeycomb  to  that  of  the  equivalent  plate  will 
therefore  be  an  indication  of  the  gain  in  strength. 

By  using  the  formulation  outlined  in  Chapter  4,  in 
conjunction  with  the  above  equations,  the  local  and  global 
deflections  of  the  structure  can  be  determined  and,  thus,  the 
corresponding  stresses  can  be  readily  calculated  by  applying 
basic  plate  theory  (a  detailed  discussion  of  the  relationship 
between  the  deflection  and  the  stress  has  been  given  by 
Timoshenko[3 ] ) . Certainly,  the  resulting  stresses  of  both 
the  top  plate  element  and  the  bottom  plate  element  are  the  sum 
of  the  stresses  corresponding  to  the  local  and  global 
deflections  (the  principle  of  superposition) . A plot  of  the 
ratio  of  the  maximum  normal  stress  of  the  top  plate  of  the 
structure  to  the  maximum  normal  stress  of  the  equivalent  plate 
is  shown  in  figure  7.6.  It  can  be  seen  from  the  plot  that 
there  is  a minimum  for  the  curve  of  the  stress  ratio, 
indicating  that  with  a specific  thickness  of  the  plate 
elements  of  the  honeycomb  structure,  there  is  a maximum  gain 
in  strength.  In  this  case,  the  optimum  stress  ratio  is 
approximately  0.5159  and  the  optimum  thickness  of  the  plate 
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Figure  7.6  Ratio  of  the  maximum  stresses 
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elements  is  approximately  0.018.  For  the  considered 
structure,  when  h=0.018,  the  maximum  normal  stress  of  the 
honeycomb  structure  is  about  52%  of  that  of  the  equivalent 
plate. 

If,  instead  of  considering  strength,  the  rigidity  is  more 
important,  then  the  main  concern  is  the  ratio  of  the  maximum 
overall  deflection  of  the  honeycomb  structure  to  that  of  a 
regular  plate  with  uniform  thickness  and  having  the  same  size. 
The  thickness  of  the  regular  plate  is  so  determined  that  its 
volume  is  the  same  as  the  honeycomb.  By  using  the  modified 
Rayleigh-Ritz  method  in  calculating  the  maximum  deflections 
of  the  honeycomb  structure  and  the  regular  plate,  the  result 
shown  in  figure  7.7  is  obtained.  The  figure  shows  a plot  of 
the  ratio  of  the  maximum  deflections  of  the  structure  and  the 
plate.  In  this  case,  an  optimum  point  can  also  be  obtained. 
This  point  corresponds  approximately  to  h=0.016,  and  the  ratio 
of  maximum  deflection  at  the  optimum  point  is  approximately 
0.192.  Accordingly,  the  maximum  deflection  of  the  optimum 
honeycomb  structure  is  19.2%  of  the  regular  plate  when  acted 
upon  by  the  same  load.  This  is  a significant  gain  in 
rigidity. 

It  seems  that  the  optimum  points  for  both  strength  and 
rigidity  are  located  close  to  each  other,  and  it  may  be 
roughly  concluded  that  for  the  particular  structure  considered 
in  this  section,  the  thickness  of  the  plate  elements  should 
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be  chosen  from  between  0.015  and  0.02  to  gain  the  most  in  both 
strength  and  rigidity. 

Since  it  is  impossible  to  have  a structure  which  is 
ideally  light  and  strong  at  the  same  time,  another  way  to 
apply  the  concept  of  optimum  design  is  to  construct  an 
objective  function  to  determine  the  optimum  thickness  of  the 
plate  elements  so  that  a compromise  between  the  rigidity  and 
the  weight  can  be  accomplished.  The  objective  function  is 
constructed  as  follows: 

F = ( 5s  / ) - Z ( V / ) (7.12) 
where  denotes  the  maximum  overall  deflection  of  the 
honeycomb  structure  (obtained  using  the  modified  Rayleigh-Ritz 
method) , denotes  the  maximum  deflection  of  the 
corresponding  solid  plate,  which  can  also  be  obtained  by  the 
proposed  method,  V is  the  total  volume  of  the  honeycomb 
structure,  is  the  volume  of  the  corresponding  solid  plate, 
and  Z is  a weighting  factor  defining  the  compromise  between 
the  volume  and  the  rigidity. 

By  using  the  modified  Rayleigh-Ritz  method  in  the 
calculation  process,  the  following  results  are  obtained. 
Figure  7.8  shows  three  curves  of  F,  corresponding  to  Z=0.5, 
1,  and  2,  plotted  versus  the  thickness  of  the  plate  elements 
of  the  honeycomb.  It  shows  that  there  is  always  an  optimum 
point  for  the  objective  function  F.  This  optimum  point 
represents  the  best  compromise  between  the  requirements  of 
less  volume  and  more  rigidity  for  each  particular  Z.  In 
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figure  7.9,  the  optimum  thickness  of  the  plate  elements  of  the 
honeycomb  is  shown  with  respect  to  different  weighting 
factors.  As  expected,  when  Z increases,  the  optimum  value  of 
h decreases.  It  can  be  seen  that  when  volume  is  important, 
the  thickness  of  the  plate  elements  of  the  honeycomb  is 
decreased  to  satisfy  the  requirement  of  less  volume.  This  is 
done,  however,  at  the  expense  of  the  rigidity  of  the 
structure.  The  maximum  overall  deflection  of  the  structure 
increases  as  Z increases  (figure  7.10).  Figure  7.11  shows  two 
curves  of  the  volume  ratio.  The  symbol  denotes  the  volume 
of  the  corresponding  solid  plate.  The  symbol  Vg  is  the  volume 
of  the  equivalent  plate  for  which  the  length  and  width  are  the 
same  as  the  honeycomb,  and  the  maximum  deflection  is  the  same 
as  the  maximum  overall  deflection  of  the  honeycomb.  It  can 
be  seen  that  it  is  advantageous  to  use  a honeycomb  structure 
to  support  transverse  loads,  since  the  volume  of  the  honeycomb 
structure  is  much  smaller  than  that  of  a solid  plate  of  the 
same  size.  It  is  also  smaller  than  that  of  an  equivalent 
plate  which  will  deflect  the  same  amount  when  acted  upon  by 
the  same  load. 

There  are  many  possibilities  for  constructing  an 
objective  function  in  accordance  with  various  design 
requirements.  The  objective  function  described  above  only 
serves  as  an  example  to  demonstrate  the  applicability  of  the 
proposed  method  to  the  analytical  aspect  of  the  optimum  design 
of  elastic  structures.  The  advantage  of  the  modified 
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Rayleigh-Ritz  method  over  the  finite  element  method  in  the 
sense  of  efficiency  can  be  seen  from  the  fact  that  for  each 
of  the  optimizations  described  above,  approximately  10 
objective  function  evaluations  were  required,  and  the 
computing  time  required  for  the  proposed  method  to  calculate 
the  deflection  of  an  element  was  observed  to  be  up  to  5 times 
faster  than  the  finite  element  method.  Furthermore,  the 
formulation  of  the  proposed  method  is  less  complicated  than 
the  finite  element  method,  and  the  computer  coding  for  the 
proposed  method  is  easier  to  achieve  than  for  the  finite 
element  method.  The  proposed  method  offers  designers  a more 
efficient  and  effective  alternative  to  handle  complex  elastic 
structures  with  various  design  requirements  and,  thus,  various 
objective  functions. 


0.600 


161 


J 


Figure  7.8  Plot  of  the  objective  function  with  respect 
to  the  thickness  of  the  plate  elements 
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Figure  7 . 10  Maximum  deflection  of  the  honeycomb  structure 
with  respect  to  different  weighting  factors 


164 


Oipj  SLUniOA 


Figure  7 . 11  Volume  ratio  with  respect  to  weighting  factors 


CHAPTER  8 

CONCLUSIONS  AND  RECOMMENDATIONS 
Conclusions 

It  can  be  concluded  from  this  study  that  the  Rayleigh- 
Ritz  method  implemented  with  undetermined  multipliers  provides 
an  excellent  tool  for  the  analysis  and  optimum  design  of  a 
wide  variety  of  elastic  structures.  The  following  is  a 
summary  of  some  of  the  findings. 

(1)  The  proposed  method  can  be  used  to  treat  almost  all 
types  of  beam  problems  provided  that  the  bending  rigidity  is 
constant.  This  includes  beams  with  different  boundary 
conditions,  beams  with  geometric  constraints,  as  well  as  beams 
acted  upon  by  different  types  of  loads.  The  main  advantage 
is  that  the  boundary  conditions  can  be  changed  without  major 
modification  of  the  formulation.  Interior  constraints  can  be 
incorporated,  and,  again,  the  formulation  remains  similar. 

(2)  The  proposed  method  can  be  applied  to  many  plate 
problems  with  different  boundary  conditions,  including  mixed 
boundary  conditions.  Plates  of  irregular  shape,  stiffened 
plates  and  honeycomb  type  structures  can  be  dealt  with  using 
the  same  formulation  as  that  used  to  handle  simple  problems. 
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(3)  The  proposed  method  can  be  applied  to  shell  problems 
for  analyzing  the  bending  effect  on  deflection.  Such  problems 
include  infinitely  long  shells  held  rigidly,  infinitely  long 
shells  reinforced  by  rings,  shells  with  fixed  ends  and  shells 
reinforced  with  both  rings  and  longitudinal  stiffeners. 

(4)  It  may  be  concluded  that  a unified  formulation  is 
possible  for  treating  a wide  variety  of  elastic  structures. 

Beside  the  concise  and  elegant  nature  of  the  formulation 
procedure,  the  solution  is  also  very  efficient  compared  to  the 
finite  element  method.  The  computer  time  reguired  to  solve 
problems  of  this  kind  by  the  finite  element  method  can  be  up 
to  several  times  that  needed  for  solution  by  the  developed 
approach.  For  the  case  of  a plate  with  an  edge  clamped  and 
the  other  three  edges  simply  supported  (as  shown  in  Chapter 
4) , the  CPU  time  for  the  proposed  method  is  96  seconds  on  an 
80286  microcomputer,  and  that  for  the  finite  element  method 
with  400  elements  is  510  seconds.  The  gain  in  efficiency  of 
the  proposed  method  over  the  finite  element  method  is  obvious. 

Another  advantage  of  the  modified  Rayleigh-Ritz  method 
is  that  in  contrast  to  the  discrete  data  provided  by  the 
finite  element  method,  the  proposed  method  generates  a 
continuous  solution. 

An  important  conclusion  of  this  study  is  that  the 
deflection  of  a plate  with  a hole  can  be  obtained  from  a plate 
without  any  hole  using  the  proposed  method  and  an  iteration 
procedure.  This  configuration  occurs  in  a variety  of 
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constructions,  and  designers  have  so  far  relied  almost 
entirely  on  the  finite  element  method  for  solutions.  In  view 
of  the  great  effort  needed  for  the  preprocessing  and 
postprocessing  of  the  finite  element  method,  the  above 
approach  can  serve  as  a valuable  tool  in  the  analysis  of 
plates  with  holes. 

In  view  of  the  advantages  of  the  proposed  method  in  the 
analysis  of  elastic  structural  elements,  the  method  can  also 
be  an  excellent  tool  for  optimum  design.  As  shown  in  the 
previous  chapter,  by  choosing  a suitable  objective  function 
and  using  the  Newton-Raphson  iteration  method,  the  optimum 
solution  can  be  obtained  in  few  iterations. 

Selection  of  the  Coordinate  Functions 

The  solution  generated  by  the  proposed  method  is  assumed 
to  be  a linear  combination  of  a set  of  predetermined 
coordinate  functions.  The  quality  of  the  result  is  highly 
dependent  upon  an  appropriate  selection  for  the  coordinate 
functions.  Some  guidance  for  selecting  the  coordinate 
functions  is  discussed  below. 

Since  the  essence  of  the  proposed  method  is  to  force 
unsatisfied  boundary  conditions  to  be  satisfied  by  choosing 
the  displacement  parameters,  it  can  be  seen  that  while  it  is 
possible  to  have  the  summation  of  the  coordinate  functions 
vanish  at  a boundary  by  adjusting  the  parameters,  it  is  not 
possible  to  make  the  summation  of  the  coordinate  functions 
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equal  to  a nonzero  value  at  a boundary  if  all  the  coordinate 
functions  vanish  at  the  boundary.  This  is  the  first  guide  in 
selecting  appropriate  coordinate  functions. 

Beside  the  above  consideration,  care  should  be  exercised 
in  the  choice  of  the  coordinate  functions  in  order  to  improve 
the  efficiency  of  the  modified  Rayleigh-Ritz  method.  For 
example,  when  sine  functions  are  chosen  as  the  coordinate 
functions,  if  the  geometry  of  the  object  is  symmetrical  about 
the  centroid,  then  choosing  odd  sine  functions  and  discarding 
even  functions  will  considerably  reduce  the  computing  time. 

Although  only  sine  and  cosine  functions  are  used  in  this 
study,  it  should  be  noted  that  other  functions,  such  as 
Chebyshev  polynomials  and  Jacobi  polynomials,  can  also  be 
taken  as  the  coordinate  functions. 

Limitations  of  the  Method 

The  modified  Rayleigh-Ritz  method,  although  possessing 
valuable  advantages  in  analyzing  and  designing  elastic 
structural  elements,  has  the  following  disadvantages: 

(1)  It  is  difficult  for  the  proposed  method  to  handle 
beams  with  variable  bending  rigidity  or  plates  with  variable 
thickness.  The  difficulty  arises  from  the  integration  of  the 
potential  energy.  It  is  difficult  to  have  a closed  form 
solution  of  the  integration  if  the  bending  or  flexure  rigidity 
is  a function  of  the  independent  variable (s). 
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(2)  Care  must  be  taken  in  selecting  the  approximation 
series  and  checking  convergency.  The  selection  of  the  type 
of  solution  function  can  be  critical  in  complex  situations. 

(3)  Applications  to  nonlinear  problems  can  be  difficult. 

(4)  The  proposed  method  is  not  readily  applicable  to 
optimum  designs  where  the  thickness  of  a plate  or  the  cross 
section  of  a beam  is  the  design  variable. 

Recommendations 

The  following  studies  are  recommended  for  future  work  so 
that  the  method  presented  here  can  be  fully  developed: 

(1)  It  is  suggested  that  a procedure  to  apply  the 
proposed  method  to  nonlinear  problems,  such  as  thick  plates 
or  plates  with  large  deflections,  be  developed. 

(2)  It  is  also  suggested  that  a procedure  to  handle 
plates  with  variable  thickness  or  beams  with  variable  bending 
rigidity  be  developed. 

(3)  It  is  further  suggested  that  a procedure  to  analyze 
and  design  a structure  which  is  an  integration  of  different 
types  of  structural  elements  be  developed. 

(4)  Finally,  as  the  original  Ritz  method  was  primarily 
developed  for  the  analysis  of  the  vibrations  of  plates[45], 
it  is  suggested  that  a procedure  for  the  dynamic  analysis  of 
a complex  structure  with  complex  boundary  conditions  could  be 
investigated,  and  a study  on  applying  the  dynamic  analysis  to 
the  design  of  elastic  structures  be  conducted. 


APPENDIX 

PROGRAM  FOR  THE  CASE  OF  HONEYCOMB  SANDWICH  PLATES 


# include 
# include 
# include 
# include 


<stdio.h> 

<math.h> 

<stdlib.h> 

<alloc.h> 


#define  PI 
#define  XSTF 
#define  YSTF 
#define  NSTF 


3.14159265358979324 

1 

3 

3 


extern  unsigned  _stklen  = 40000; 


short 

short 

short 

short 

short 

short 

short 

short 

double 

double 

double 

double 

double 

double 

double 

double 

double 

double 

double 

double 

double 

double 

double 

double 

double 

main() 

{ 

short 

short 


TENF,  CENF,  BENF,  PENF,  NENF,  NINS=2 ; 

XSEC=9 ; 

YSEC=5 ; 

NODE  ; 

TERM=15 ; 

TERM2 ; 

*npes ; 

*inst ; 

epss=l . Oe-12 ; 

volume,  toln; 

qt=1.0,  qb=0.00; 

lx=3.0,  ly=1.0,  tp; 

ts=0. 010000,  h=0. 100000,  hi; 

Ep=3.0e07,  Gp=11.54e06,  mup=0.3,  density=0 . 283 
yield_stress=6 . 0e04 , ultimate_stress=l . 0e05 ; 
dpr,  dpc,  El,  GJ; 

Inst [NSTF] ; 

EIv[NSTF] ; 
ornt[NSTF] ; 
xstr[NSTF] ; 
ystr[NSTF] ; 

*u , *v ; 

*eta ; 

*lgmr ; 

*dm; 


loopl,  loop2 ; 
indexl ; 
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double  *cdpm,  *tdpm,  *bdpm; 

void  ini_sf(); 
void  ini_uv() ; 
void  ini_eta() ; 
void  ini_din(); 
void  ini_EIv() ; 
void  solve_LM() ; 
void  cmpt_dfpm() ; 
void  cmpt_defl ( ) ; 
void  memo_free ( ) ; 
short  *s_mem_alc() ; 
double  *meino_allc  ( ) ; 

scant ( "%le” , &tp  ); 

TERM2  = TERM  * TERM; 

printf ( ” Basic  information. \n"  ); 

printf(  '•  1.  MATERIAL:  \n''  ); 

printf ( " Poison's  ratio:  %.3f\n",  mup  ); 

printf ( " Yang's  modulus:  %.3e\n",  Ep  ) ; 

printf ( "\n  2.  LOAD:\n"  ); 

printf ( " load  on  top:  %.2f,  load  on  bottom:  %.2f\n", 

qt,  qb  ) ; 


npes  = s_mem_alc(  NSTF,  "npes"  ) ; 
inst  = s_mem_alc(  2*NINS,  "inst"  ) ; 
ini_sf(  ornt,  Inst,  xstr,  ystr,  npes  ); 

TENF  = 0; 

for(  loopl=0;  loopKNSTF;  ++loopl  ) 

TENF  +=  ( ( *(npes+loopl)  + 1 ) / 2 ); 

BENF  = TENF; 

CENF  = TENF; 

PENF  = TENF  + CENF  + BENF; 

NENF  = PENF  + NINS ; 

eta  = memo_allc(  NENF,  "eta"  ) ; 
ini_eta ( eta  ) ; 

u = memo_allc(  NENF,  "u"  ) ; 

V = memo_allc(  NENF,  "v"  ) ; 

Igmr  = memo_allc(  NENF,  "Igmr"  ) ; 

toln  = (double) XSTF  * lx  + (double) YSTF  * ly; 
hi  = h - (double) 2.0  * tp; 
volume  = (double) 2.0  * tp  * lx  * ly; 
volume  +=  toln  * hi  * ts; 

El  = Ep  * ts  * pow(  hi,  (double)!  ) / (double) 12 . 0 ; 
ini_EIv(  EIv  ) ; 

dpr  = Ep  * pow(  tp,  (double)!  ) 

/ ( (double) 1.0  - mup  * mup  ); 
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dpr  /=  (double) 12 . 0 ; 

dpc  = ( hi  + tp  ) / (double) 2 . 0 ; 

dpc  = Ep  * tp  * pow(  dpc,  (double) 2 ) ; 

dpc  /=  ( (double) 1.0  - mup  * mup  ) ; 

dpc  +=  dpr; 


printf(  "\n  3.  PLATE;\n''  ); 

printf(  ” length:  %.2f,  width:  %.2f,  thickness  %f\n”, 

lx,  ly,  tp  ) ; 

printf(  ” distance  between  plates  %.2f\n",  hi  ); 
printf ( " flextural  rigidity  (reference)  %e\n",  dpr  ); 

printf ( ” flextural  rigidity  (nutral  pi)  %e\n",  dpc  ) ; 

printf ( ”\n  4.  STIFFENER”  ) ; 

printf ( ” (based  on  a quarter  of  the  plate) \n”  ) ; 

printf ( ” no  of  stiffener:  %d\n” , NSTF  ) ; 

printf ( ” reference  bending  rigidity;  %e\n”,  El  ); 

printf ( ” thickness  of  stiffener:  %f\n",  ts  ); 

printf ( " Tp  / Ts  = %e\n” , tp  / ts  ) ; 

printf ( ” In  on  xs  ys  rd  cp\n”  ) ; 

for(  loopl=0;  loopKNSTF;  ++loopl  ) 

printf(  ” #%2d  %.2f  %.2f  %.2f  %.2f  %.2f  %2d\n”  , 

loopl,  * (Inst+loopl) , * (ornt+loopl) , * (xstr+loopl) , 

* (ystr+loopl) , * (EIv+loopl) , * (npes+loopl)  ); 


ini_uv(  u,  V,  inst,  ornt,  xstr,  ystr,  Inst,  eta  ) ; 


printf ( ”\n  locations  of  reinforced  points\n"  ); 

printf ( ” no  x y eta\n"  ) ; 

for(  loopl=0;  loopKNENF;  ++loopl  ) 
printf ( ” %d  %f  %f  %f\n”, 

loopl,  *(u+loopl),  *(v+loopl),  *(eta+loopl)  ); 


printf ( ”\n  5.  INTERSECTION: \n”  ); 
for(  loopl=0;  loopKNINS;  ++loopl  ) 

printf ( ” #%d:  stf  #%d  and  stf  #%d\n", 

loopl,  * (inst+loopl*2) , * (inst+loopl*2+l)  ); 


dm  = memo_allc(  TERM2 , "dm”  ); 
ini_dm(  dm,  TERM  ) ; 
cdpm  = memo_allc(  TERM2,  "cdpm”  ); 
tdpm  = memo_allc(  TERM2,  "tdpm”  ); 
bdpm  = memo_allc(  TERM2 , "bdpm”  ); 

solve_LM(  Igmr,  TERM  ) ; 

cmpt_dfpm(  cdpm,  tdpm,  bdpm,  TERM  ) ; 

cmpt_def 1 ( cdpm,  tdpm,  bdpm,  TERM  ); 

memo_free(  bdpm,  "bdpm”  ); 
memo_free(  tdpm,  "tdpm”  ); 
memo_free(  cdpm,  "cdpm”  ) ; 


) 
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void  ini_sf(  ornt,  Inst,  xstr,  ystr,  npes  ) 
double  *ornt,  *lnst; 

double  *xstr,  *ystr; 

short  *npes ; 

{ 

short  loopl ; 

ornt[0]  = (double) 0 . 0 ; 
ornt[l]  = PI  / (double) 2.0; 
ornt [2]  = PI  / (double) 2 . 0; 
lnst[0]  = lx; 
lnst[l]  = ly; 

Inst [2]  = ly; 
xstr[0]  = (double) 0.0; 
ystr[0]  = ly  / (double) 2.0; 
xstr[l]  = lx  / (double) 4.0; 
ystr[l]  = (double) 0.0; 
xstr [2]  = lx  / (double) 2.0; 
ystr [2]  = (double) 0.0; 

for(  loopl=0;  loopK (XSTF+l)/2 ; ++loopl  ) 
npes [loopl]  = 15; 

for(  loopl=(XSTF+l)/2 ; loopK (XSTF+YSTF+2 ) /2 ; ++loopl  ) 
npes [loopl]  = 10; 


} 


void  ini_xy(  x,  y ) 
double  *x,  *y; 

{ 

short  loop; 
double  ds; 

ds  = lx  / lx  / (double) 2.0  / (double) ( XSEC  - 1 ); 
loop  = 0 ; 

do  *(x+loop)  = (double) loop  * ds; 
while ( ++loop  < XSEC  ) ; 

ds  = ly  / ly  / (double) 2.0  / (double) ( YSEC  - 1 ); 
loop  = 0; 

do  *(y+loop)  = (double) loop  * ds; 
while ( ++loop  < YSEC  ) ; 

} 


void  ini_eta(  eta  ) 
double  *eta; 

{ 

short  loopl,  loop2; 


short  indexl=0,  index2=0; 
double  deta; 

for(  loopl=0;  loopKTENF;  ++loopl  ) 

*(eta+loopl)  = (double) 0.0; 

for(  loopl=0;  loopKNSTF;  ++loopl  ) 

{ 

* (npes+loopl)  +=  1; 

deta  = * (Inst+loopl)  / (double) ( * (npes+loopl)  ); 

* (npes+loopl)  /=  2; 
index2  +=  * (npes+loopl) ; 

for(  loop2=indexl;  loop2<index2 ; ++loop2  ) 

*(eta+loop2)  = (double) ( loop2  - indexl  + 1 ) * 

indexl  = index2 ; 

} 


} 


void  ini_uv(  u,  v,  inst,  ornt,  xstr,  ystr,  Inst,  eta  ) 
short  *inst; 
double  *u,  *v; 

double  *ornt,  *xstr,  *ystr,  *lnst; 
double  *eta ; 

{ 

short  loopl,  loop2 ; 
short  indexl=0,  index2=0; 
double  dcos,  dsin,  dtan; 
double  cos(),  sin(),  tan(),  sqrt()  ; 

for(  loopl=0;  loopKNSTF;  ++loopl  ) 

{ 

dcos  = cos(  * (ornt+loopl)  ); 
dsin  = sin(  * (ornt+loopl)  ); 

if  ( fabs(  dcos  ) < epss  ) dcos  = (double) 0.0; 
if  ( fabs(  dsin  ) < epss  ) dsin  = (double) 0.0; 
index2  +=  * (npes+loopl) ; 

for(  loop2=indexl ; loop2<index2 ; ++loop2  ) 

{ 

*(u+loop2)  = * (xstr+loopl)  + *(eta+loop2)  * dcos; 
*(u+loop2)  /=  lx; 

* (u+loop2+CENF)  = *(u+loop2); 

* (u+loop2+CENF+TENF)  = *(u+loop2); 

*(v+loop2)  = * (ystr+loopl)  + *(eta+loop2)  * dsin; 
*(v+loop2)  /=  ly; 

* (v+loop2+CENF)  = *(v+loop2); 

* (v+loop2+CENF+TENF)  = *(v+loop2); 

*(eta+loop2)  /=  * (Inst+loopl) ; 

* (eta+loop2+CENF)  = * (eta+loop2 ) ; 
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deta ; 
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* (eta+loop2+CENF+TENF)  = * (eta+loop2 ) ; 

} 

indexl  = index2 ; 

} 

indexl  = PENF; 
index2  = 0; 

for(  loopl=0;  loopKNSTF;  ++loopl  ) 

{ 

for(  loop2=loopl+l ; loop2<NSTF;  ++  loop2  ) 

{ 

if  ( fabs(  * (ornt+loopl)  - * (ornt+loop2)  ) > l.Oe-6  ) 

{ 

if  ( fabs(  fabs(  * (ornt+loopl)  ) 

- PI  / (double) 2.0  ) <=  l.Oe-07  ) 
dcos  = * (xstr+loopl) ; 

else 

{ 

if  ( fabs(  fabs(  * (ornt+loop2)  ) 

- PI  / (double) 2.0  ) <=  l.Oe-07  ) 
dcos  = * (xstr+loop2) ; 

else 

{ 

dcos  = * (ystr+loop2)  - * (ystr+loopl) ; 
dcos  +=  ( tan(  * (ornt+loopl)  ) * * (xstr+loopl)  ); 
dcos  -=  ( * (xstr+loop2)  * tan(  * (ornt+loop2 ) ) ); 
dcos  /=  ( tan(  * (ornt+loopl)  ) 

- tan(  * (ornt+loop2)  ) ); 

} 

} 

if  ( dcos  <=  lx  / (double) 2.0  &&  dcos  >=  0.0  ) 

{ 

dsin  = tan(  * (ornt+loopl)  ) 

* ( dcos  - * (xstr+loopl)  ) ; 

dsin  +=  * (ystr+loopl) ; 

if  ( dsin  <=  ly  / (double) 2.0  &&  dsin  >=  0.0  ) 

{ 

*(u+indexl)  = dcos  - * (xstr+loopl) ; 

*(v+indexl)  = dsin  - * (ystr+loopl) ; 

*(u+indexl)  *=  *(u+indexl); 

*(v+indexl)  *=  *(v+indexl); 

* (eta+indexl)  = 

sqrt(  *(u+indexl)  + *(v+indexl)  ); 
*(u+indexl)  = dcos  - * (xstr+loop2 ) ; 

*(v+indexl)  = dsin  - * (ystr+loop2) ; 

*(u+indexl)  *=  *(u+indexl); 

*(v+indexl)  *=  *(v+indexl); 

*(v+indexl)  = sqrt(  *(u+indexl)  + *(v+indexl)  ); 
*(v+indexl)  /=  * (lnst+loop2) ; 

*(u+indexl)  = * (eta+indexl)  / * (Inst+loopl) ; 
indexl++ ; 

* (inst+index2)  = loopl; 


176 


index2++; 

* (inst+index2)  = loop2 ; 
index2++; 

) 

} 

} 

} 

} 


) 


void  ini_EIv(  EIv  ) 
double  *EIv; 

{ 

short  loopl ; 
double  epss=l . Oe-08 ; 

for(  loopl=0;  loopKNSTF;  ++loopl  ) 
if  ( fabs(  * (ornt+loopl)  ) < epss 

&&  fabs(  * (ystr+loopl)  - ly  / (double) 2.0  ) <=  epss  ) 
*(EIv+loopl)  = El  / (double) 2.0; 
else 

if  ( fabs(  * (ornt+loopl)  - PI  / (double) 2.0  ) < epss 
&&  fabs(  * (xstr+loopl)  - lx  / (double) 2.0  ) <=  epss  ) 
*(EIv+loopl)  = El  / (double) 2.0; 
else 

*(EIv+loopl)  = El; 


) 


void  ini_dm(  dm,  terml  ) 
double  *dm; 
short  terml ; 

{ 

double  tmp; 

short  loopl , loop2 ; 

short  indexl , index2 ; 

for(  loopl=0;  loopKterml;  ++loopl  ) 

{ 

tmp  = ( (double) ( 2 * loopl  + 1 ) / lx  ) ; 
tmp  *=  tmp; 

indexl  = loopl  * terml; 

for(  loop2=0;  loop2<terml;  ++loop2  ) 

{ 

index2  = indexl  + loop2 ; 

*(dm+index2)  = ( (double) ( 2 * loop2  + 1 ) / ly  ) ; 
*(dm+index2)  *=  * (dm+index2) ; 
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* (dm+index2)  +=  tmp; 

* (dm+index2)  *=  * (dm+index2 ) ; 

} 


} 


double  qmn(  m,  n,  i ) 
short  m,  n,  i; 

{ 

double  tmpl , tmp2 ; 
double  sin(); 

tmpl  = (double) ( 2 * m + 1 ) * PI; 

tmp2  = (double) ( 2 * n + 1 ) * PI; 

tmpl  = sin(  tmpl  * *(u+i)  ); 

if  ( fabs(  tmpl  ) < epss  ) tmpl  = 

tmp2  = sin(  tmp2  * *(v+i)  ); 

if  ( fabs(  tmp2  ) < epss  ) tmp2  = 

return ( tmpl  * tmp2  ) ; 

} 


double  rmn(  k,  i ) 
short  k,  i; 

{ 

double  tmpl; 
double  sin() ; 

tmpl  = (double) ( 2 * k + 1 ) * PI; 
tmpl  = sin(  tmpl  * *(eta+i)  ); 
if  ( fabs(  tmpl  ) < epss  ) tmpl  = 

return ( tmpl  ) ; 

} 


void  solve_LM(  Igmr,  TERM  ) 
double  *lgmr; 
short  TERM ; 

{ 

double  *a; 
short  ele  a; 


(double) 0.0; 
(double) 0.0; 


(double) 0.0; 


double  *memo_allc ( ) ; 
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void  iiiemo_free()  ; 
void  cmpt_a() ; 
void  cmpt_lginr  ( ) ; 
short  Choleski ( ) ; 

ele_a  = NENF  * NENF  + NENF; 
ele_a  /=  2 ; 

a = memo_allc(  ele_a,  "a”  ) ; 

cmpt_a(  a,  ele_a,  TERM  ) ; 

cmpt_lginr(  Igmr,  TERM  ) ; 

if  ( ! Choleski ( a,  Igmr,  NENF  ) ) 

{ 

printf  ("Fail  to  perform  Choleski 's  decomposition  (1)  .\n")  > 
printf ( "Program  terminates. \n"  ); 
exit(O) ; 

} 

memo_free(  a,  "a"  ) ; 

) 


void  cmpt_a(  a,  ele_a,  TERM  ) 

double  *a; 

short  ele_a,  TERM; 

{ 

double  cons,  cont,  tmpO,  tmpl,  tmp2,  tmp3 ; 
short  idnl , idn2 ; 

short  loopO,  loopl,  loop2,  loop3 , loop4 ; 
short  indexl,  index2,  index3,  index4,  indexS; 
short  indexs , indexe ; 
double  gmn() , rmn() ; 

for(  loopl=0;  loopl<ele_a;  ++loopl  ) 

*(a+loopl)  = (double) 0 . 0 ; 
index3  = 0 ; 

for(  loopl=l;  loopK=CENF;  ++loopl  ) 
index3  +=  loopl; 
index4  = 0 ; 

for(  loopl=l;  loopK=CENF+TENF;  ++loopl  ) 
index4  +=  loopl; 

cont  = (double) 1.0  / lx  / ly; 
cons  = cont  / ( dpc  - dpr  ) ; 
cont  *=  ( (double) 2.0  / dpr  ); 
indexl  =0; 

for(  loopl=0;  loopKCENF;  ++loopl  ) 

{ 
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index!  +=  CENF; 

index4  +=  ( CENF  + TENF  ) ; 

indexl  +=  loopl; 

for(  loop2=0;  loop2<=loopl ; ++loop2  ) 

{ 

if  ( ( index2=indexl+loop2  ) > ele_a  ) 

{ 

printf 

("Intend  to  use  memory  outside  the  allocated  area.\n") ; 
printf ( "Program  terminates. \n"  ) ; 
exit(O) ; 

} 

*(a+index2)  = (double) 0.0; 
tmpO  = (double) 0.0; 
loop!  = 0; 
do 

{ 

indexS  = loop!  * TERM; 
loop4  = 0; 
do 

{ 

tmpl  = qmn(  loop!,  loop4,  loopl  ); 
tmpl  *=  qmn(  loop!,  loop4,  loop2  ); 
tmpl  /=  * (dm+index5+loop4) ; 
tmpO  +=  tmpl; 

) 

while ( ++loop4  < TERM  ) ; 

} 

while ( ++loop!  < TERM  ) ; 

*(a+index2)  +=  ( tmpO  * cons  ); 
tmpl  = tmpO  * cont; 

* (a+index!+index2)  = *(a+index2)  + tmpl; 

* (a+index!+index2-CENF)  = *(a+index2); 

* (a+index4+index2)  = *(a+index2)  + tmpl; 

} 

} 

index2  = 0 ; 

for(  loopl=CENF;  loopKCENF+TENF ; ++loopl  ) 

{ 

indexl  +=  loopl; 

index!  = indexl  + index2 ; 

for(  loop2=loopl-CENF+l ; loop2<CENF;  ++loop2  ) 

{ 

index!  +=  ( CENF  + loop2  ) ; 

* (a+indexl+loop2)  = *(a+index!); 

} 

++index2 ; 

} 

index!  = 0 ; 

for(  loopl=CENF+l;  loopK=CENF+TENF;  ++loopl  ) 
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index3  +=  loopl; 

for(  loopl=CENF+TENF;  loopKPENF;  ++loopl  ) 

{ 

indexl  +=  loopl; 

for(  loop2=0;  loop2<CENF;  ++loop2  ) 

{ 

index2  = indexl  + loop2 ; 

*(a+index2)  = * (a+index2-index3) ; 

* (a+index2+CENF)  = *(a+index2); 

) 

index3  +=  CENF; 

} 

indexl  = 0; 
index3  = 0 ; 
index4  = 0 ; 

for(  loop4=0;  loop4<NSTF;  ++loop4  ) 

{ 

cons  = pow(  * (lnst+loop4) , (double)!  ) / * (EIv+loop4) ; 
cons  /=  (double) 2.0; 
index4  +=  * (npes+loop4) ; 

for(  loopl=index3 ; loopl<index4;  ++loopl  ) 

{ 

indexl  +=  loopl; 

for(  loop2=index3 ; loop2<=loopl ; ++loop2  ) 

{ 

if  ( ( index2=indexl+loop2  ) > ele_a  ) 

{ 

printf 

("Intend  to  use  memory  outside  the  allocated  area.\n") ; 
printf ( "Program  terminates . \n"  ); 
exit(O) ; 

} 

tmpO  = (double) 0.0; 
loop!  = 0; 
do 

{ 

tmpl  = rmn(  loop!,  loopl  ); 
tmpl  *=  rmn(  loop!,  loop2  ); 
tmpl  /= 

pow(  (double) ( 2 * loop!  + 1 ),  (double) 4 ); 
tmpO  +=  tmpl; 

} 

while ( ++loop!  < TERM  ) ; 

*(a+index2)  +=  ( tmpO  * cons  ); 

} 

} 

index!  = index4 ; 

} 


indexl 


0; 
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for(  loopl=0;  loopKPENF;  ++loopl  ) 
indexl  +=  loopl; 

for(  loopl=PENF;  loopKPENF+NINS ; ++loopl  ) 

{ 

indexl  +=  loopl; 
loop3  = loopl  - PENF; 
idnl  = * (inst+2*loop3) ; 
idn2  = * (inst+2*loop3+l) ; 
indexs  = 0; 

for(  loop3=0;  loop3<idnl;  ++loop3  ) 
indexs  +=  * (npes+loop3) ; 
indexe  = indexs  + * (npes+idnl) ; 

cons  = pow(  *(lnst+idnl)  ,(double)3  ) / *(EIv+idnl); 
cons  /=  (double) 2.0; 

for(  loop2=indexs;  loop2<=loopl  &&  loop2<indexe; 

++loop2  ) 

{ 

if  ( ( index2=indexl+loop2  ) > ele_a  ) 

{ 

printf 

("Intend  to  use  memory  outside  the  allocated  area.\n"  ) ; 
printf ( "Program  terminates. \n"  ); 
exit(O) ; 

} 

tmpO  = (double) 0.0; 
loop4  = 0; 
do 

{ 

tmp3  = (double) ( 2 * loop4  + 1 ) * PI ; 
tmpl  = sin(  tmp3  * *(u+loopl)  ); 
if  ( fabs(  tmpl  ) < epss  ) tmpl  = (double) 0.0; 
tmpl  *=  rmn(  loop4,  loop2  ); 

tmp2  = pow(  (double) ( 2 * loop4  + 1 ),  (double) 4 ); 

tmpO  +=  ( tmpl  / tmp2  ) ; 

} 

while ( ++loop4  < TERM  ) ; 
tmpO  *=  cons; 

*(a+index2)  +=  tmpO; 

} 

for(  loop3=idnl;  loop3<idn2;  ++loop3  ) 
indexs  +=  * (npes+loop3) ; 
indexe  = indexs  + * (npes+idn2) ; 

cons  = pow(  *(lnst+idn2)  , (double)!  ) / *(EIv+idn2); 
cons  /=  (double) 2.0; 
cons  = -cons; 

for(  loop2=indexs;  loop2<=loopl  &&  loop2<indexe; 

++loop2  ) 

{ 

if  ( ( index2=indexl+loop2  ) > ele_a  ) 
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{ 

printf 

("Intend  to  use  memory  outside  the  allocated  area.\r^"  ) ' 
printf ( "Program  terminates. \n"  ); 
exit(O) ; 

} 

tmpO  = (double) 0.0; 
loop4  =0; 
do 

{ 

tmp3  = (double) ( 2 * loop4  + 1 ) * PI ; 
tmpl  = sin(  tmp3  * *(v+loopl)  ) ; 
if  ( fabs(  tmpl  ) < epss  ) tmpl  = (double) 0.0; 
tmpl  *=  rmn(  loop4,  loop2  ); 

tmp2  = pow(  (double) ( 2*  loop4  + 1 ),  (double) 4 ); 
tmpO  +=  ( tmpl  / tmp2  ) ; 

} 

while ( ++loop4  < TERM  ) ; 
tmpO  *=  cons; 

*(a+index2)  +=  tmpO; 

} 

for(  loop2=PENF;  loop2<=loopl ; ++loop2  ) 

{ 

loop3  = loop2  - PENF; 

if  ( ( index2=indexl+loop2  ) > ele_a  ) 

{ 

printf 

("Intend  to  use  memory  outside  the  allocated  area.\n") ; 
printf ( "Program  terminates. \n"  ); 
exit(O) ; 

} 

if  ( * (inst+loop3*2)  ==  idnl  ) 

{ 

cons  = pow(  * (Inst+idnl) , (double)!  ) / *(EIv+idnl); 

cons  /=  (double) 2.0; 

tmpO  = (double) 0.0; 

loop4  = 0; 

do 

{ 

tmp3  = (double) ( 2 * loop4  + 1 ) * PI ; 

tmpl  = sin(  tmp3  * *(u+loopl)  ); 

if  ( fabs(  tmpl  ) < epss  ) tmpl  = (double) 0.0; 

tmp2  = sin(  tmp3  * *(u+loop2)  ); 

if  ( fabs(  tmp2  ) < epss  ) tmp2  = (double) 0.0; 

tmpl  *=  tmp2 ; 

tmp2  = pow(  (double) ( 2 * loop4  + 1 ) , (double) 4 ) ; 
tmpO  +=  ( tmpl  / tmp2  ) ; 

) 

while ( ++loop4  < TERM  ) ; 

*(a+index2)  +=  tmpO  * cons; 

} 
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if  ( * (inst+loop3*2+l)  ==  idn2  ) 

{ 

cons  = pow(  * (lnst+idn2) , (double) 3 ) / *(EIv+idn2); 

cons  /=  (double) 2.0; 

tmpO  = (double) 0.0; 

loop4  = 0 ; 

do 

{ 

tmp3  = (double) ( 2 * loop4  + 1 ) * PI ; 

tmpl  = sin(  tmp3  * *(v+loopl)  ); 

if  ( fabs(  tmpl  ) < epss  ) tmpl  = (double) 0.0; 

tmp2  = sin(  tmp3  * *(v+loop2)  ); 

if  ( fabs(  tmp2  ) < epss  ) tmp2  = (double) 0.0; 

tmpl  *=  tmp2 ; 

tmp2  = pow(  (double)  ( 2 * loop4  + 1 ) , (double) 4 ) ; 
tmpO  +=  ( tmpl  / tmp2  ) ; 

} 

while ( ++loop4  < TERM  ) ; 

*(a+index2)  +=  tmpO  * cons; 

} 

} 

} 

) 


void  cmpt_lgmr(  Igmr,  TERM  ) 
double  *lgmr; 
short  TERM ; 

{ 

short  loopl,  loop2,  loop3 , loop4 ; 
short  indexl , index2 , index! ; 
double  tmpl,  tmp2 , tmp3 , tmp4 ; 
double  qmn() ; 

loopl  = 0; 

while ( loopl  < NENF  ) * (lgmr+ (loopl++) ) = (double) 0.0; 

index2  = CENF; 
index!  = CENF; 

for(  loop4=0;  loop4<NSTF;  ++loop4  ) 

{ 

index!  +=  * (npes+loop4) ; 

for(  loop3=index2 ; loop3<index3 ; ++loop3  ) 

{ 

for(  loopl=0;  loopKTERM;  loopl++  ) 

{ 

tmpl  = (double) ( 2 * loopl  + 1 ) * PI; 

indexl  = loopl  * TERM; 

for(  loop2=0;  loop2<TERM;  loop2++  ) 

{ 

tmp2  = (double) ( 2 * loop2  + 1 ) * PI ; 


tmp2  *=  tmpl; 

tmp4  = qmn(  loopl,  loop2,  loop3  ); 

* (lgmr+loop3)  += 

( tinp4  / tmp2  / * (dm+indexl+loop2 ) ) 

) 

} 

* (lgmr+loop3+TENF)  = * (lgmr+loop3) ; 

} 

index2  = index3 ; 

} 

tmpl  = (double) 2.0  * qt  / dpr; 
tmp2  = (double) 2.0  * qb  / dpr; 
loop3  = 0; 
do 

{ 

tmp3  = (double) 0.0; 

if  ( loop3  >=  CENF  ) tmp3  = tmpl; 

if  ( loop3  >=  (TENF+CENF)  ) tmp3  = tmp2 ; 

if  ( loop3  >=  PENF  ) tmp3  = (double) 0.0; 

* (lgmr+loop3)  *=  tmp3 ; 

} 

while ( ++loop3  < NENF  ) ; 

} 


short  Choleski(  A,  lambda,  max_lambda  ) 
double  *A,  * lambda; 
short  max  lambda ; 


{ 

void 

forward_substitute ( ) ; 

void 

backward_substitute ( ) ; 

short 

i=0,  j=0,  k,  1,  m,  n; 

do 

{ 

1 = ( j +=  i/  j+i  ) ; 

for  ( k=0;  k<i;  ++k  ) 

*(A+1)  -=  pow(  *(A+j+k),  (double)2  ); 
if  ( *(A+1)  <=  0.0  ) 

{ 

printfC'Stop  %d  %le\n",  1,  *(A+1)  ); 
return (0) ; 

} 

*(A+1)  = sqrt(  *(A+1)  ); 
n = 1; 

for  ( k=i+l;  k<max_lambda ; ++k  ) 
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{ 

n +=  k; 

for  ( m=0;  m<i;  ++in  ) 

*(A+n)  -=  *(A+n-m-l)  * * (A+l-m-1) ; 
*(A+n)  /=  *(A+1); 

} 

} while  ( ++i  < max_laitibda  ) ; 

forward_substitute ( A,  lambda,  max_lambda  ) ; 
backward_substitute ( A,  lambda,  max_lambda  ) ; 

return ( 1 ) ; 

} 


void  forward_substitute ( A,  lambda,  max_lambda  ) 
double  *A,  *lambda; 
short  max_lambda; 

{ 

short  i , j , k ; 
k = 0 ; 

for  ( i=0;  i<max_lambda ; ++i  ) 

{ 

k +=  i; 

for  ( j=0;  j<i;  ++j  ) 

*(lambda+i)  -=  *(lambda+j)  * *(A+k+j); 
*(lambda+i)  /=  *(A+k+i); 

} 


} 


void  backward_substitute ( A,  lambda,  max_lambda  ) 
double  *A,  * lambda; 
short  max_lambda; 

{ 

short  i=0,  j=0,  k,  1,  m,  n; 


— max_lambda ; 

for  ( i=max_lambda;  i>=0;  — i ) 

{ 

k = 0; 

for  ( 1=1;  l<=i;  ++1  ) k +=  1; 
m = 0; 

for  ( j=i+l;  j<=max_lambda;  ++j  ) 

{ 

m +=  j ; 

*(lambda+i)  -=  *(lambda+j)  * *(A+k+i+m); 

} 
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* (lambda+i)  /-  *(A+k+i); 

} 

} 


double  *memo_allc ( size,  string  ) 
char  *string; 

short  size; 

{ 

double  *p; 

void  err_meiiio  ( ) ; 

size  *=  sizeof (double) ; 

p = (double  *)malloc(  (unsigned) size  ); 
if  ( p ==  NULL  ) err_memo(  string  ) ; 

return ( p ) ; 

} 


short  *s_mem_alc(  size,  string  ) 
char  * string; 

short  size; 

{ 

short  *p ; 

void  err_memo() ; 

size  *=  sizeof (short) ; 

p = (short  *)iiialloc(  (unsigned)  size  ); 
if  ( p ==  NULL  ) err_meino(  string  ); 

return ( p ) ; 

) 


void  meino_free(  array_po inter , string  ) 
char  * string; 

double  *array_pointer ; 

{ 

free(  array_pointer  ) ; 

printf ( "Array  %s  is  released\n" , string  ) ; 


} 
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void  err_memo(  string  ) 
char  *string; 

{ 


printf ("Warning! !\n") ; 

printf("Fail  to  allocate  memory  for  array 

printf ( "%s",  string  ); 

print f ( "\nProgram  terminates . \n" ) ; 

exit (0) ; 


} 


void  cmpt_dfpm(  cdpm,  tdpm,  bdpm,  TERM  ) 
double  *cdpm,  *tdpm,  *bdpm; 
short  TERM ; 

{ 

double  tmpl,  tmp2 , tmp3 , tmp4 ; 
short  loopl,  loop2 , loop3 ; 
short  indexl , index2 ; 
double  qmn() ; 

tmp3  = lx  * ly  / PI  / PI; 
tmpl  = tmp3  * qt; 
tmp3  *=  qb; 

for(  loopl=0;  loopKTERM;  loopl++  ) 

{ 

indexl  = loopl  * TERM; 

tmp2  = tmpl  / (double) ( 2 * loopl  + 1 ) 

tmp4  = tmp3  / (double) ( 2 * loopl  + 1 ) 

for(  loop2=0;  loop2<TERM;  loop2++  ) 

{ 

index2  = indexl  + loop2 ; 

* (cdpm+index2)  = (double) 0.0; 

* (tdpm+index2)  = tmp2  / (double) ( 2 * 

* (bdpm+index2)  = tmp4  / (double) ( 2 * 

} 

} 

for(  loopl=0;  loopKTERM;  ++loopl  ) 

{ 

indexl  = loopl  * TERM; 

for(  loop2=0;  loop2<TERM;  ++loop2  ) 

{ 

index2  = indexl  + loop2 ; 
tmpl  = (double) 0.0; 
loop3  = 0 ; 
do 

{ 

tmp2  = qmn(  loopl,  loop2,  loop3  ); 
tmpl  +=  ( * (lgmr+loop3)  * tmp2  ); 


loop2  + 1 ) ; 
loop2  + 1 ) ; 
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} 

while ( ++loop3  < CENF  ); 

* (cdpm+index2)  +=  tmpl; 

tmpl  = (double) 0.0; 
loop3  = CENF; 
do 

{ 

tmp2  = qmn(  loopl,  loop2,  loop3  ); 
tmpl  +=  ( * (lgmr+loop3)  * tmp2  ); 

} 

while ( ++loop3  < TENF+CENF  ); 

* (tdpm+index2)  -=  tmpl; 

* (cdpm+index2)  +=  tmpl; 

tmpl  = (double) 0.0; 
loop 3 = CENF+TENF; 
do 

{ 

tmp2  = qmn(  loopl,  loop2,  loop3  ); 
tmpl  +=  ( * (lgmr+loop3)  * tmp2  ); 

) 

while ( ++loop3  < PENF  ) ; 

* (bdpm+index2)  -=  tmpl; 

* ( cdpm+index2 ) +=  tmpl ; 

) 


tmpl  = (double) 8.0  / lx  / ly  / ( PI  * PI  * PI  * PI  ) ; 
tmp2  = tmpl  / ( dpc  - dpr  ) ; 
tmpl  *=  ( (double) 2.0  / dpr  ); 
for(  loopl=0;  loopKTERM;  ++loopl  ) 

{ 

indexl  = loopl  * TERM; 

for(  loop2=0;  loop2<TERM;  ++loop2  ) 

{ 

index2  = indexl  + loop2 ; 
tmp3  = tmpl  / * (dm+index2) ; 

* (cdpm+index2)  *=  ( tmp2  / *(dm+index2)  ); 

* (tdpm+index2)  *=  tmp3 ; 

* (bdpm+index2)  *=  tmp3 ; 

) 


} 


void  cmpt_defl(  cdpm,  tdpm,  bdpm,  TERM  ) 
double  *cdpm,  *tdpm,  *bdpm; 
short  TERM; 

{ 
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FILE  *fpgt,  *fpot,  *fpgb,  *fpob,  *fpgc, 
double  *cdef,  *tdef,  *bdef; 
double  *x,  *y; 

double  tmpl,  tmp2,  tmp3 , tmp4 ; 
short  loopl,  loop2,  loop3,  loop4 ; 
short  indexl,  index2,  index3 ; 
void  ini_xy() ; 
void  outp_grd() ; 
void  memo_free() ; 
double  *memo_allc() ; 
double  sin() ; 

NODE  = XSEC  * YSEC; 

X = meino_allc(  XSEC,  "x"  ); 
y = meino_allc(  YSEC,  "y”  ) ; 
ini_xy(  x,  y ) ; 

cdef  = memo_allc(  NODE,  "cdef”  ) ; 
tdef  = memo_allc(  NODE,  "tdef”  ) ; 
bdef  = menio_allc ( NODE,  "bdef"  ) ; 

for(  loopl=0;  loopKNODE;  ++loopl  ) 

{ 

* (cdef+loopl)  = (double) 0 . 0 ; 

* (tdef+loopl)  = (double) 0.0; 

* (bdef+loopl)  = (double) 0 . 0 ; 

) 

for(  loopl=l;  loopKXSEC;  ++loopl  ) 

{ 

indexl  = loopl  * YSEC; 

for(  loop2=l;  loop2<YSEC;  ++loop2  ) 

{ 

loop3  = 0; 
do 

{ 

index2  = loop3  * TERM; 

tmp4  = (double) ( 2 * loop3  + 1 ) * 

tmpl  = sin(  tmp4  * *(x+loopl)  ); 

if  ( fabs(  tmpl  ) < epss  ) tmpl  = 

loop4  = 0 ; 

do 

{ 

tmp4  = (double) ( 2 * loop4  + 1 ) 
tmp4  = sin(  tmp4  * *(y+loop2)  ); 
if  ( fabs(  tmp4  ) < epss  ) tmp4 
tmp4  *=  tmpl ; 
index3  = index2  + loop4 ; 
tmp2  = * (cdpm+ index! ) * tmp4 ; 
tmp3  = * (tdpm+ index!)  * tmp4 ; 
tmp4  *=  * (bdpm+index3) ; 
index!  = indexl  + loop2 ; 


*fpoc; 


PI; 

(double) 0.0; 

* PI; 

= (double) 0.0; 
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* (cdef+index3)  +=  tmp2; 

* (tdef+index3)  +=  tmp3; 

* (bdef+index3)  +=  tmp4 ; 

} 

while ( ++loop4  < TERM  ) ; 

} 

while ( ++loop3  < TERM  ) ; 

} 

} 

fpot  = fopen(  "eye2t.out" , "w”  ) 
fpgt  = fopen(  "eye2t.grd” , ”w”  ) 
fpob  = fopen(  "eye2b.out",  "w"  ) 
fpgb  = fopen(  ”eye2b.grd” , ”w"  ) 
fpoc  = fopen(  "eye2c. out” , ”w”  ) 
fpgc  = fopen(  ”eye2c.grd”,  ”w”  ) 
outp_grd(  fpgt,  fpot,  tdef  ) ; 

outp_grd(  fpgb,  fpob,  bdef  ) ; 

outp_grd(  fpgc,  fpoc,  cdef  ) ; 

f close ( fpgc  ) ; 
f close ( fpoc  ) ; 
fclose(  fpgb  ) ; 
f close ( fpob  ) ; 
f close ( fpgt  ) ; 
fclose(  fpot  ) ; 

memo_free(  cdef,  "cdef”  ) ; 
memo_free(  bdef,  "bdef”  ) ; 
memo_free(  tdef,  "tdef”  ) ; 
memo_free(  y,  ”y”  ); 
itiemo_free  ( x,  ”x"  ) ; 

} 


void  outp_grd(  fpg,  fpo,  df  ) 

FILE  *fpg/  *fpo; 

double  *df; 

{ 

short  loopl,  loop2 ; 
short  indexl , index2 ; 

for(  loopl=0;  loopKXSEC;  ++loopl  ) 

{ 

indexl  = loopl  * YSEC; 

for(  loop2=0;  loop2<YSEC;  ++loop2  ) 

{ 

index2  = indexl  + loop2 ; 

fprintf(  fpo,  ” w[%2d] [%2d]=%le\n” , 

loopl,  loop2,  *(df+index2)  ); 

} 

} 
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loopl  = 2 * XSEC  - 1; 
loop2  = 2 * YSEC  - 1; 

fprintf(  fpg,  "DSAA\n"  ); 
fprintf(  fpg,  "%d  %d\n",  loopl,  loop2  ); 
fprintf(  fpg,  ”0.0  %.2f\n”,  lx  ); 

fprintf(  fpg,  ”0.0  %.2f\n”,  ly  ); 

fprintf(  fpg,  ”0.0  1.0\n”  ); 

indexl  = YSEC  - 1; 

index2  = NODE  - YSEC; 

for(  loopl=0;  loopK indexl;  ++loopl  ) 

{ 

for(  loop2=loopl;  loop2<index2 ; loop2+=YSEC  ) 
fprintf(  fpg,  ”%le\n”,  *(df+loop2)  ); 
for(  ; loop2>=loopl ; loop2-=YSEC  ) 

fprintf ( fpg,  ”%le\n”,  *(df+loop2)  ); 

} 

for(  ; loopl>=0;  — loopl  ) 

{ 

for(  loop2=loopl;  loop2<index2 ; loop2+=YSEC  ) 
fprintf ( fpg,  "%le\n”,  *(df+loop2)  ); 
for(  ; loop2>=loopl;  loop2-=YSEC  ) 

fprintf ( fpg,  ”%le\n”,  *(df+loop2)  ); 

} 


} 
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